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Preface 


The  thesis  is  the  second  in  a  series  AFIT  thesis  research  efforts  to  develop  a  new 
post-processing  Kalman  filter  for  the  Completely  Integrated  Reference  Instrumentation 
System  (CIRIS).  The  ultimate  goal  is  the  “Advanced  CIRIS  Filter,”  designed  around 
the  current  CIRIS  II  LN-39  INS,  which  processes  both  standard  and  differential  mode 
Global  Positioning  System  (GPS)  measurements  in  addition  to  the  current  ground-based 
transponder  measurements.  The  new  Kalman  filter  is  required  to  increase  the  position  and 
velocity  estimation  accuracy  of  CIRIS  so  that  CIRIS  will  remain  a  more  accurate  estimator 
relative  to  other  types  of  navigation  systems.  This  research  is  sponsored  by  the  Central 
Inertial  Guidance  Test  Facility  (CIGTF),  6585th  Test  Group,  Holloman  AFB,  New  Mexico. 

The  new  CIRIS  Kalman  filter  developed  in  this  and  previous  research  is  built  upon  the 
Multimode  Simulation  for  Optimal  Filter  Evaluation  (MSOFE)  software  developed  by  the 
Avionics  Laboratory  at  Wright-Patterson  AFB,  Ohio.  Although  MSOFE  was  conceived 
and  implemented  as  a  simulation  tool,  its  thoughtfully  designed  structure  makes  it  readily 
adaptable  for  use  with  real  measurements.  It  seems  to  be  a  consensus  among  the  AFIT 
students  who  have  used  MSOFE  that  it  is  a  difficult  tool  to  master.  I  agree,  but  having 
become  a  journeyman  MSOFE  user  myself,  I  see  that  the  long  learning  curve  is  rooted  in 
the  difficult  concepts  of  Kalman  filter  theory.  MSOFE  itself  is  very  flexible  and  obviously 
well  engineered.  In  an  academic  setting,  MSOFE  is  exactly  the  type  of  tool  needed  in  a 
“Kalman  filter  laboratory”  to  make  the  concept  of  a  Kalman  filter  come  to  life.  Hopefully, 
this  thesis  shows  that  MSOFE  can  also  be  put  to  good  use  outside  of  research  laboratories 
and  academic  institutions. 

At  the  onset  of  my  particular  research  I  had  hoped  to  pick  up  at  the  point  where 
my  predecessor,  Capt  Joseph  Solomon,  had  left  off,  and  to  concentrate  on  developing  GPS 
error  models  and  extending  the  filter  structure  to  incorporate  real  GPS  measurements. 
However,  the  idea  of  “serial”  research,  especially  in  a  graduate-school  setting,  is  often  over- 
optimistic.  I  had  to  acquire  and  assimilate  a  great  deal  of  totally  new  (to  me)  knowledge 
before  I  could  comprehend  what  Capt  Solomon  had  accomplished  and  what  remained  to 
be  done.  This  need  to  “get  up  to  speed,”  combined  with  the  late  availability  of  real  GPS 
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data  and  my  personal  interest  in  the  discipline  of  software  engineering,  led  me  to  focus 
my  efforts  on  small  refinements  to  the  error  model  on  which  the  filter  is  based,  an  almost 
complete  rewrite  of  the  filter  software,  and  additional  analysis  of  the  filter’s  performance 
using  regular  CIRIS  (non-GPS)  data  sets.  I  hope  that  my  forthcoming  assignment  to 
CIGTF  will  permit  me  to  continue  the  development  of  the  Advanced  CIRIS  Filter  in 
parallel  with  follow-on  research  at  AFIT. 

I  am  indebted  to  several  individuals  for  their  assistance  during  my  stay  at  AFIT.  My 
thesis  advisor,  Lt  Col  Zdzislaw  (Stan)  Lewantowicz,  guided  me  through  this  effort  with 
much  patience  and  good  advice.  If  any  individual  can  be  said  to  have  a  natural  aptitude 
for  the  mathematics  of  error  model  development,  it  is  Col  Lewantowicz.  I  was  privileged  to 
learn  the  theory  of  stochastic  estimation  from  Dr.  Peter  Maybeck.  Dr.  Maybeck's  ability  as 
an  instructor  and  the  depth  of  his  knowledge  both  exceed,  in  stochastics  terminology,  the 
+3<r  threshold.  At  the  Avionics  Lab,  Mr.  Stanton  Musick,  the  prime  force  behind  MSOFE, 
and  his  assistant,  Mr.  Robert  Urbanic,  provided  exceptional  assistance  during  the  time  I 
was  learning  to  use  MSOFE.  At  CIGTF,  Mr.  Gordon  Simkin  and  Mr.  Francisco  Ramerez 
answered  my  many  questions  and  provided  me  with  CERIS  data.  I  also  wish  to  thank 
Capt  Joseph  Solomon,  Dr.  Robert  Ewing,  and  Mr.  Don  Smith.  Last,  because  his  name 
starts  with  a  “Z,”  but  first  in  my  book,  is  Mr.  Daniel  Zambon,  director  of  AFIT’s  Signal 
Information  Processing  Lab.  Mr.  Zambon  is  genuinely  dedicated  to  providing  reliable 
computer  resources  in  support  of  student  research  efforts.  He  patiently  assisted  me  on 
many  occasions  and  for  this  I  am  extremely  appreciative. 

I  also  wish  to  acknowledge  the  lifelong  support  of  my  parents,  Charles  Edward 
Snodgrass  and  Shirley  Ann  Snodgrass.  Their  love  and  encouragement  is  priceless. 


Faron  Britt  Snodgrass 
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Abstract 

/ 

The  Completely  Integrated  Reference  Instrumentation  System  (CIRIS)  is  operated 
by  the  Central  Inertial  Guidance  Test  Facility  (CIGTF),  at  Holloman  AFB,  New  Mex¬ 
ico.  CIRIS  functions  as  a  reference  navigation  system  used  for  evaluating  the  accuracy 
of  other  types  of  navigations  systems.  As  a  reference  standard,  the  root-mean-square 
errors  in  CIRIS  estimates  of  aircraft  trajectory  variables  must  be  maintained  an  order  of 
magnitude  smaller  than  the  corresponding  errors  of  the  navigation  system  under  test.  The 
primary  hardware  components  of  CIRIS  are  a  reference  inertial  navigation  system  (INS), 
a  baro-altimeter,  an  array  of  ground-based  transponders,  a  transponder  interrogator,  and 
data  recording  equipment.  The  transponder  equipment  provides  transponder-to-aircraft 
range  and  range-rate  measurements  during  test  flights.  The  primary  software  compo¬ 
nent  of  CIRIS  is  a  Kalman  filter  program  which  processes  the  recorded  measurements  and 
estimates  the  true  position  and  velocity  of  the  aircraft  throughout  the  test  flight.  ^ 

The  estimation  accuracy  of  CIRIS  must  be  increased  so  that  CIRIS  can  serve  as  a 
benchmark  for  measuring  the  accuracy  of  Global  Positioning  System  (GPS)  aided  inertial 
navigation  systems.  This  thesis  documents  the  continuation  of  research  to  develop  com¬ 
pletely  new  Kalman  filter  software  for  CIRIS.  A  70-state  filter,  based  on  the  Multimode 
Simulation  for  Optimal  Filter  Evaluation  (MSOFE)  program,  developed  in  previous  re¬ 
search  is  the  starting  point.  This  70-state  filter  models  error  dynamics  associated  with  the 
CIRIS  Litton  LN-39  INS,  baro-altimeter,  and  transponder  equipment.  In  this  thesis,  the 
model  for  atmospheric  effects  on  transponder  range  measurements  is  refined  and  the  filter 
is  modified  to  process  barometric  altitude  measurements  in  addition  to  the  transponder 
measurements.  The  performance  of  the  resulting  filter  is  evaluated  using  real  CIRIS  data 
recorded  during  a  slow  speed  ground  test  and  an  aircraft  flight  test.  The  filter  position  and 
velocity  estimates  are  compared  to  independent  measurements  of  the  same  quantities.  The 
structure  for  a  companion  fixed-interval  smoother  program  is  proposed  and  discussed,  but 
not  implemented.  Future  research  is  expected  to  extend  the  filter  to  process  differential 
mode  GPS  measurements. 
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Continued  Development  and  Analysis  of  a  New  Extended 
Kalman  Filter  for  the  Completely  Integrated  Reference 
Instrumentation  System  (CIRIS) 


I.  Introduction 

This  thesis  describes  the  continuing  development  and  analysis  of  an  extended  Kalman 
filter  for  the  Advanced  Completely  Integrated  Reference  Instrumentation  System  (CIRIS). 
The  major  goals  of  this  thesis  effort  are  filter  software  development,  filter  model  valida¬ 
tion  via  performance  evaluation,  and  description  of  a  “smoothing”  algorithm  that  makes 
maximum  use  of  the  filter  calculations. 

CIRIS  is  a  system  for  evaluating  the  accuracy  of  added  or  unaided  inertial  navigation 
systems  (INS).  It  is  operated  by  the  Central  Inertial  Guidance  Test  Facility  (CIGTF), 
6585th  Test  Group,  Air  Force  Systems  Command,  Holloman  AFB,  New  Mexico.  Its  pri¬ 
mary  hardware  components  are  a  Litton  LN-39  INS,  an  array  of  ground  based  transpon¬ 
ders,  and  a  transponder  interrogator  [2,11].  The  transponders  are  located  at  precisely 
surveyed  points  throughout  the  continental  United  States.  For  a  CIRIS  test  flight,  the 
INS  being  evaluated  (test  article)  is  placed  in  an  aircraft  alongside  the  LN-39  INS  and  the 
transponder  interrogator.  During  the  flight  of  the  aircraft,  the  interrogator  requests  and 
receives  range  and  velocity  information  from  selected  transponders.  The  transponders  are 
interrogated  one  at  a  time.  Typically,  a  “window”  containing  eight  or  fewer  transponders 
is  selected  and  the  interrogator  repeatedly  cycles  through  only  these  transponders.  The 
specific  transponders  contained  in  the  window  may  change  during  the  flight.  The  range, 
velocity,  baro-altimeter,  and  LN-39  INS  information  is  currently  processed  by  an  extended 
Kalman  filter  algorithm  running,  in  real-time,  on  a  digital  computer.  The  Kalman  filter 
provides  accurate  estimates  of  the  aircraft  position,  velocity,  and  attitude  throughout  the 
aircraft  flight.  This  information  is  used  as  a  reference  against  which  the  position,  velocity, 
and  attitude  calculations  of  the  test  article  are  compared. 
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1.1  Background 

The  original  CIRIS  system  was  developed  by  CIGTF  in  the  early  1970’s  and  became 
operational  in  1975  [4].  It  uses  a  Litton  LN-15  INS  which  is  no  longer  routinely  supported 
by  Air  Force  maintenance  depots.  This  original  system,  called  CIRIS  I,  is  still  in  use.  A 
newer  version  of  CIRIS,  based  on  the  LN-39  INS,  is  called  CIRIS  II.  Both  versions  of  CIRIS 
are  based  on  a  14-state  extended  Kalman  filter  of  which  only  11  states  are  actually  used. 
States  1-9  are  the  primary  INS  position,  velocity,  and  attitude  error  states.  States  10  and 
11  model  errors  in  the  baro- altimeter  aiding  of  the  vertical  channel.  States  12-14  were 
originally  intended  to  model  errors  in  Doppler  velocity  aiding  but  are  not  used  because 
initial  tests  indicated  the  Doppler  Eliding  did  not  provide  a  significant  increase  in  accuracy. 
All  subsequent  references  to  CIRIS  refer  to  CIRIS  II.  The  term  “AdvEinced  CIRIS”  is  used 
here  to  refer  to  the  high-order  (50  or  more  states)  post-processing  filters  developed  in  this, 
previous,  and  subsequent  AFIT  thesis  research. 

The  11-state  CIRIS  Kalman  filter  runs  on  a  Hewlett-Packard  HP  1000  computer 
carried  aboard  the  aircraft.  The  position  and  velocity  estimates  are  recomputed  every 
second.  The  current  CIRIS  I  system  provides  latitude  and  longitude  accurate  to  13  feet 
(ft)  l<r,  altitude  accurate  to  40  ft  1<t,  north  and  west  velocity  accurate  to  0.1  feet/second 
( fps )  1  <r,  and  vertical  velocity  accurate  to  0.4  fps  ler  [4].  The  l<r  standcird  deviation  value 
is  for  a  Gaussian  (normal)  error  distribution. 

Currently,  the  position  and  velocity  estimates  provided  by  CIRIS  are  more  accurate 
than  those  from  almost  any  other  type  of  navigation  system.  Thus  CIRIS  is  a  reference 
standard  for  evaluating  the  accuracies  of  navigation  systems.  However,  more  accurate 
navigation  systems  are  becoming  available  due  to  increased  manufacturing  precision,  new 
sensor  technologies,  more  powerful  navigation  computers,  Eind  satellite  based  navigation 
systems.  Since  it  is  essential  that  the  reference  navigation  system  be  at  least  an  order  of 
magnitude  more  accurate  the  navigation  systems  being  evaluated,  the  accuracy  of  CIRIS 
must  be  increased. 

A  powerful  characteristic  of  a  properly  designed  Kalman  filter  is  the  ability  to  process 
measurements  of  some  quantity  from  two  or  more  independent  sources  and  produce  an 
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estimate  of  that  quantity  that  is  more  accurate  than  the  individual  measurements.  This 
characteristic  depends  on  the  accuracy  of  the  mathematical  models  of  the  physical  system 
as  the  basis  for  proper  Kalman  filter  design. 

Three  ways  of  increasing  CIRIS  accuracy  are  apparent.  The  first  approach  is  to  ex¬ 
pand  the  number  of  states  in  the  CIRIS  Kalman  filter  to  model  additional  time-correlated 
error  sources  in  the  existing  LN-39  INS  and  transponder  measurement  equipment.  This 
is  the  approach  taken  by  Captain  Joseph  Solomon  in  his  AFIT  thesis  research  [13].  A 
second,  complementary  approach  is  to  reprocess  the  filter  estimates  with  a  type  of  back¬ 
ward  running  optimal  estimator  known  as  a  smoother.  The  third  approach  is  to  aid  CIRIS 
with  Global  Position  System  (GPS)  measurements.  The  continuation  of  the  first  approach 
and  the  initial  investigation  of  the  second  approach  are  the  primary  goals  of  this  thesis 
research.  The  third  approach  will  be  the  subject  of  future  AFIT  thesis  research. 

1.2  Research  Objectives 

There  are  three  main  research  objectives.  The  first  is  to  improve  the  filter  software 
developed  by  Solomon.  This  includes  improving  the  source  code  structure  and  efficiency, 
making  the  filter  easier  to  use,  and  adding  features  and  capabilities. 

The  second  objective  is  to  evaluate,  using  real  CIRIS  data,  the  performance  of  a 
refined  version  of  the  70-state  filter  developed  by  Solomon.  This  includes  limited  tuning 
of  selected  filter  parameters  to  improve  performance. 

The  third  objective  is  the  development  of  a  fixed-interval  smoothing  algorithm.  This 
type  of  smoother  is  an  optimal  estimator  that  processes  previously  filtered  information 
backwards  in  time.  This  requires  that  certain  data  from  the  forward  pass  be  recorded  for 
subsequent  use  by  the  smoother.  The  smoothing  algorithm  starts  with  the  most  recent 
data  and  proceeds  backward  to  the  beginning  of  the  interval.  This  characteristic  prevents 
the  use  of  the  smoother  in  real-time;  it  must  be  used  as  a  data  post-processor.  Since  the 
CIGTF  mission  does  not  normally  require  read-time  data  processing,  the  operation  of  the 
advanced  CIRIS  filter  amd  the  smoother  as  post-processors  will  enhance  the  quality  of  the 
reference  trajectory. 
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1.3  Research  Approach 

The  performance  of  the  revised  70-state  filter  algorithm  is  determined,  and  the  error 
model  validated,  using  data  from  a  series  of  tests  conducted  at  the  6585^  Test  Group’s 
rocket  sled  test  track.  An  instrument  rack  containing  the  CIRIS  LN-39  INS  and  transpon¬ 
der  interrogator  was  installed  in  a  Recovery  Support  Vehicle  (RSV).  The  RSV  was  towed 
along  the  track  at  speeds  of  40-60  fps  while  the  CIRIS  data  was  recorded.  At  the  same 
time,  the  RSV’s  position  along  the  track  was  recorded  by  the  independent  track  data  acqui¬ 
sition  system.  Both  sets  of  data  were  time-tagged  using  an  Inter-Range  Instrumentation 
Group  (IRIG)  time  standard  common  to  both  data  acquisition  systems.  The  relatively 
slow  speeds  are  chosen  to  make  the  data  collection  time  interval  reasonably  large.  Since 
the  starting  location  is  known  precisely  and  the  relative  locations  of  the  track  position  sen¬ 
sors  are  known  within  a  millimeter  [7],  the  track  measurements  are  an  accurate  standard 
by  which  the  advanced  CIRIS  filter  performance  can  be  evaluated.  Filter  model  validation 
is  achieved  by  analysis  of  measurement  residual  characteristics  and  by  comparison  of  the 
filter’s  position  find  velocity  estimates  to  those  derived  from  the  track  measurements. 

The  accuracy  of  the  70-state  filter  estimates  is  also  evaluated  using  CIRIS  data 
from  an  actual  aircraft  flight.  The  reference  standard  for  this  test  is  the  measurements  of 
aircraft  position  and  velocity  provided  by  the  laser  tracking  system  at  the  Yuma,  Arizona, 
test  range.  The  average  accuracy  of  these  laser  measurements  is  roughly  twice  that  of 
CIRIS  [7]. The  CIRIS  II  filter  estimates  are  also  compared  to  the  laser  measurements  so 
that  accuracy  of  Advanced  CIRIS  relative  to  CIRIS  II  may  be  evaluated. 

The  smoothing  algorithm  suggested  by  this  thesis  research  is  a  discrete-time  inverse 
covariance  formulation.  The  algorithm  was  originally  developed  by  Meditch  [10]  and  fur¬ 
ther  described  by  Maybeck  [9],  The  method  makes  use  of  the  forward  pass  calculations 
of  the  covariance  matrix  and  the  state  dynamics  matrix  for  the  backward  pass  calculation 
of  a  smoother  gain  matrix.  The  forward  pass  time-history  for  the  dynamics  matrix,  the 
covariance  matrix,  and  the  state  vector  must  be  stored  for  later  use  during  the  smoothing 
procedure. 
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1-4  Thesis  Overview 

Chapter  II  describes  the  fundamental  concepts  used  in  this  research.  This  includes 
coordinate  frame  definitions  and  transformations,  basic  Kalman  filter  theory,  shaping  fil¬ 
ters  for  time-correlated  noises,  and  measurement  residual  characteristics.  The  computer 
software  packages  used  to  support  this  research  are  described  here. 

Chapter  III  provides  a  description  of  the  error  states  and  structure  of  the  70-state 
model  on  which  the  Advanced  CERIS  filter  is  based.  The  method  of  compensating  for 
“lever  arm”  effects  is  presented  in  detail. 

Chapter  IV  describes  the  test  procedures  and  70-state  filter  performance  for  CIRIS 
data  collected  during  the  test  track  runs.  The  existence  of  am  amomalous  periodic  variation 
in  some  of  the  transponder  range  measurement  data  is  documented  here. 

Chapter  V  describes  the  70-state  filter  performamce  for  CIRIS  data  collected  during 
a  long  aircraft  flight  over  the  Yuma  test  range.  The  Advamced  CIRIS  filter  output  is 
compaired  to  position  amd  velocity  measurements  from  the  laser  tracking  equipment.  The 
accuracy  of  Advanced  CIRIS  relative  to  CIRIS  II  is  discussed  here. 

Chapter  VI  describes  the  adgorithm  for  the  proposed  fixer-interval  smoother  program 
and  evailuates  its  practicadity  for  use  with  the  MSOFE-based  Advanced  CIRIS  filter.  The 
proposed  smoothing  adgorithm  is  am  adaptation  of  a  procedure  originadly  designed  for  use 
with  filters  based  on  discrete-time  dynaimics  models. 

Chapter  VII  summaurizes  the  finad  results  of  this  research.  Conclusions  baised  on  these 
results  and  recommendations  for  future  research  in  this  airea  Eire  described  here. 
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II.  Theory 


This  chapter  describes  the  fundamental  concepts  and  theory  on  which  this  research  is 
based.  Because  of  the  commonality  between  this  research  and  that  of  Solomon,  significant 
portions  of  this  chapter  are  paraphrased  from  Solomon’s  thesis  [13]. 

2. 1  Reference  Frames 

In  terrestial  navigation  applications,  position  is  usually  specified  in  terms  of  a  spher¬ 
ical  coordinate  system  of  geodetic  latitude,  longitude,  and  altitude.  Velocity  is  usually 
specified  in  terms  of  a  rectangular  coordinate  system  based  on  two  compass  directions 
and  a  vertical  direction  such  as  north,  west,  and  up.  It  is  often  necessary  to  relate  such 
geographically  descriptive  coordinates  to  an  earth-centered,  earth-fixed  (ECEF)  rectan¬ 
gular  coordinate  system,  (E- frame).  Several  geodetic  systems  have  been  developed,  but 
the  most  accurate  is  the  World  Geodetic  System  1984  (WGS  84).  This  system  is  provided 
to  the  Department  of  Defense  (DOD)  by  the  Defense  Mapping  Agency  (DMA).  In  the 
WGS  84  system,  the  earth  is  modeled  as  an  oblate  ellipsoid  as  illustrated  in  Figure  2.1. 
The  parameters  shown  in  Figure  2.1,  the  ellipse  eccentricity  e,  and  the  ellipse  flattening 
constant  (ellipticity)  /  are  defined  in  Table  2.1  [3J.  The  Xe  axis  is  parallel  to  the  meridian 
plane  which  includes  the  Bureau  International  de  L’Heure  (BIH).  This  is  the  Zero  Meridian 
(Greenwich  Meridian).  The  Ye  axis  is  rotated  90  degrees  to  the  east  along  the  equatorial 
plane.  The  Ze  axis  is  colinear  with  the  earth’s  spin  axis  and  completes  the  right-handed, 
earth-fixed,  orthogonal  coordinate  system.  The  nonlinear  equations  [15]  relating  these 
rectangular  coordinates  to  the  geodetic  latitude  (£),  longitude  (A),  and  altitude  (h)  are 


Rn 

A 

\/l  -  e2  sin2  L 

(2.1) 

xe 

=  (i2„  +  h)  cos  L  cos  A 

(2.2) 

n 

=  (Rn  +  h)  cos  L  sin  A 

(2.3) 

=  ( Rn(  1  -  e2 )  +  h)  sin  L 

(2.4) 

where  Rn  is  the  radius  of  curvature  in  the  prime  vertical. 
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Figure  2.1.  WGS  84  Ellipsoidal  Model 


Table  2.1.  WGS  84  Parameters 


Parameter 

Definition 

Value 

X,  Y',Z< 

ECEF  Coordinate  Frame  Axes 

not  applicable 

uie 

Anguish  Rate  of  the  Earth 

E2SEZB9EBEBH1 

A 

Semimajor  Axis 
(Equatorial  Radius) 

6378137  m 

B 

Semiminor  Axis 
(Polar  Radius) 

6356752.3142  m 

e 

First  Eccentricity 

0.0818191908426 

f 

Flattening  (Ellipticity) 

0.00335281066474 

9o 

Equatorial  Acceleration 
of  Gravity 

9.7803267714  m/s2 
(32.087686258  ft/s2) 
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Figure  2.2.  Coordinate  Frame  Relationships 


The  LN-39  INS  implements  a  wander-azimuth,  local-level  platform  mechanization 
[14].  The  true  platform  frame  (T-frame)  is  illustrated  in  Figure  2.2.  The  right-handed, 
orthogonal  set  of  axes  pointing  north,  west,  and  up  are  designated  N,  W ,  and  V  respec¬ 
tively.  This  frame  is  defined  as  the  NWU  navigation  frame  (iV-frame).  The  right-handed, 
orthogonal  T-frame  axes  axe  designated  X,  Y,  and  Z  respectively;  this  is  also  referred  to 
as  the  wander-azimuth  frame.  The  wander  angle  alpha  (a)  is  the  angle  between  the  X  and 
N  axes  and  between  the  Y  and  W  axes  created  by  a  counter-clockwise  (positive)  rotation 
about  the  Z,U  axis,  of  the  T-frame  with  respect  to  the  N -frame.  When  the  wander  angle 
is  zero,  the  T-frame  axes  axe  parallel  with  the  corresponding  A-frame  axes.  In  Figure  2.2, 
geodetic  latitude,  longitude,  and  altitude  above  the  WGS  84  reference  ellipsoid  are  desig¬ 
nated  by  the  variables  L,  A,  and  h,  respectively.  These  three  parameters  define  the  vehicle 
position  in  the  navigation  frame. 

Another  useful  coordinate  frame  is  the  aircraft  body  (5-frame).  This  is  reference 
frame  that  is  “attached"  to  the  body  of  the  aircraft  so  that  the  orientation  of  the  5-frame 


2-3 


with  respect  to  the  rigid  aircraft  body  is  constant.  The  origin  of  the  5- frame  is  usually 
taken  to  be  the  center  of  the  inertial  measurement  unit  (IMU),  a  subsystem  on  the  INS.  In 
this  thesis,  the  orientation  of  the  5-frame  axes  is  defined  so  that  the  5-frame  is  parallel 
to  the  NWU  navigation  frame  when  the  aircraft  nose  is  pointed  north  with  wings  level. 
The  body  frame  Xb  axis  points  out  the  nose  of  the  aircraft,  the  Yb  axis  points  out  the  left 
wing,  find  the  Zb  axis  points  up.  This  is  not  the  usual  convention  for  a  body  frame  but  it 
is  more  convenient  in  this  application.  The  orientation  of  the  5-frame  with  respect  to  the 
N -frame  is  defined  by  the  three  Euler  angles:  heading,  pitch,  and  roll.  A  positive  heading 
angle  t/>  is  created  by  clockwise  rotation  of  the  body  about  the  Zb  axis,  a  positive  pitch 
angle  9  is  created  by  clockwise  rotation  of  the  body  about  the  Yb  axis,  and  a  positive  roll 
angle  <fr  is  created  by  a  clockwise  rotation  of  the  body  about  the  Xb  axis. 

2. 2  Coordinate  Transformations 

This  section  defines  the  coordinate  transformation  matrices  needed  for  relating  one 
reference  frame  to  another.  In  addition  to  the  E,  N,  and  T  frames,  the  LN-39  Systems 
Engineering  Analysis  Report  [14]  defines  a  platform  frame  (5-frame)  and  a  computation 
frame  (C-frame).  The  true  frame,  platform  frame,  and  computation  frame  refer  to  the 
same  general  coordinate  frame.  However  the  platform  and  computation  frames  are  slightly 
misaligned,  by  error  angles,  with  respect  to  the  true  frame  and  with  respect  to  each  other. 
The  platform  frame  is  the  wander-azimuth  orthogonal  frame  “attached”  to  the  INS  plat¬ 
form  such  that  5-frame  coincides  with  the  T-frame  only  when  the  platform  is  truly  locally 
level.  The  computation  frame  is  the  wander-azimuth  orthogonal  frame  that  is  locally  level 
at  the  latitude  and  longitude  indicated  by  the  INS  computer.  Thus  the  C-frame  coincides 
with  the  T-frame  only  if  the  INS  indicated  position  is  free  of  error.  The  transformations 
from  the  true  frame  to  the  platform  and  computation  frames  are  defined  [14]  as 

[X]P  =  [I  +  0]  [X]r  (2.5) 

0  <pz  ~<t>Y 

4>  —  -<f>z  0  4>x 

4>y  ~<fx  0 
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[X]c 


6® 


=  [I  +  60][X\T 

0  6®z  —6®y 

=  -6®z  0  S®x 

6&y  ~6®x  0 


(2.6) 


The  S ®  and  ^variables  are  skew-symmetric  matrices  representing  small  misalignment 
angles.  This  concept  is  used  by  Litton  in  the  derivation  of  the  INS  error  equations.  The 
coordinate  transformation  matrices  for  the  E,  N ,  and  P  frames  are  now  defined 


[*f 

CNP 


cNp  [xf 


(2.7) 


cos  a  -  sin  a  0 
sin  a  cos  a  0 

0  0  1 

C%[X]N 

-  sin  L  cos  A  sin  A  cos  L  cos  A 

-  sin  L  sin  A  -  cos  A  cos  L  sin  A 

cos  L  0  sin  L 

f~>E 


(2.8) 


(2.9) 


.  (2.10) 


(2.11) 


—  sin  L  cos  A  cos  a  +  sin  A  sin  a  sin  L  cos  A  sin  a  4-  sin  A  cos  a 

—  sin  L  sin  A  cos  a  —  cos  A  sin  a  sin  L  sin  A  sin  a  —  cos  A  cos  a 

cos  L  cos  a  —  cos  L  sin  a 


cos  L cos  A 
cos  L  sin  A 
sin  L 


One  more  coordinate  transformation  matrix,  Cp,  is  required  to  project  vectors  co- 
ordinatized  in  the  5-frame  onto  the  IV-frame.  That  is 


[Xf  =  C%  [Af 


(2.12) 


Because  of  the  relative  complexity  of  Cg ,  it  is  developed  here  in  some  detail.  Consider 
an  aircraft  initially  pointed  north  with  wings  level.  The  5-frame  is  then  aligned  with 
the  N -frame.  Any  arbitrary  orientation  of  the  5-frame  relative  to  the  N -frame  may  be 
described  by  an  ordered  set  of  rotations  about  the  5-frame  axes.  The  5-frame  is  first 
rotated  about  Zf,  by  the  heading  angle  ip.  This  intermediate  orientation  of  the  5-frame 
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Figure  2.3.  Body  to  Navigation  Frame  Transformation 

defines  intermediate  frame  U.  Next,  the  B-frame  is  rotated  about  Yj,  by  the  pitch  angle 
9.  This  intermediate  orientation  of  the  B-frame  defines  intermediate  frame  V .  Last,  the 
B-frame  is  rotated  about  Xi,  by  the  roll  single  <p.  This  results  in  the  new  orientation  of 
the  B-frame  relative  to  the  N -frame.  This  sequence  of  rotations  is  shown  in  Figure  2.3. 
In  terms  of  these  three  rotations,  Equation  (2.12)  may  be  expressed  as 

[X ]N  =  eg  C$  Cl  [X]B  (2.13) 

where 

Cg  =  Cy  Cy  Cg  (2.14) 

Expressed  in  terms  of  the  roll,  pitch,  and  heading  angles,  the  three  component  matrices 
are 

cos  ip  sin  ip  0 

C$  =  -  sin  ip  cos  tp  0 

0  0  1 
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2.3  Kalman  Filter  Theory 

A  Kalman  filter  is  a  recursive  data  processing  algorithm  used  to  compute  time- 
dependent  estimates  of  quantities  of  interest.  When  certain  assumptions  are  made,  the 
Kalman  filter  algorithm  has  been  mathematically  proven  to  be  optimal  with  respect  to 
several  criteria  [8].  The  term  “Kalman  filter”  is  used  to  describe  several  variations  of  a 
common  algorithmic  basis.  The  model  of  a  continuous-time  linear  system  with  discrete- 
time  linear  measurements  is  used  here  to  present  the  Kalman  filter  concept. 

The  dynamics  of  the  continuous-time  system  are  assumed  well  modeled  by  a  set  of 
n  coupled  first-order  linear  differential  equations  during  some  time  period  of  interest,  T. 
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Such  a  description  of  a  physical  system  is  termed  a  state-variable  model.  These  equations 
are  written  using  vector/matrix  notation  as 

x(t)  =  F  (t)x(t)  -f  B(t)u(0  +  G(t)w(t)  (2.18) 

where  F  is  a  n-by-n  matrix  describing  the  undriven  (homogeneous)  state  dynamics,  x  is  a 
n-by-1  state  vector,  B  is  a  n-by-r  control  distribution  matrix,  u  is  a  r-by-1  deterministic 
control  input  vector,  G  is  a  n-by-s  driving  noise  distribution  matrix,  and  w  is  a  3-by-l 
stochastic  driving  noise  vector.  The  matrices  F,  B,  G,  and  the  vector  u  are,  in  general, 
piecewise  continuous  functions  of  time.  The  elements  of  the  vector  w  are  required  to  be 
white  with  zero-mean  Gaussian  statistics.  The  term  “white”  means  uncorrela^ed  in  time. 
This  concept  is  described  mathematically  by 

E{vr{t)}  =  0  (2.19) 

£{w(t)wr(t')}  =  Q  (t)tf(t-t')  (2.20) 

where  the  notation  E{ •}  denotes  expected  value.  Matrix  Q  is  the  s-by-s  driving  noise 

« 

strength  matrix,  and  6(t  -  t')  in  the  Dirac  delta  (unit  impulse)  function.  Matrix  Q  is 
symmetric,  positive  semidefinite,  and  whose  diagonal  elements  represent  the  power  spectral 
densities,  constant  at  all  frequencies,  of  the  corresponding  scalar  elements  of  w. 

At  discrete  times  t,  £  T,  m  noise-corrupted  measurements  are  available.  These 
measurements,  assumed  to  be  linear  combinations  of  the  states  and  white  discrete-time 
measurement  noises,  are  described  by 

Zi  =  z  (<j)  =  H(<,-)x(<i)  +  \(U)  (2.21) 

where  z  is  the  m-by-1  measurement  vector,  H  is  the  m-by-n  measurement  (state  obser¬ 
vation)  matrix,  x  is  the  n-by-1  state  vector,  and  v  is  the  m-by-1  stochastic  measurement 
noise  vector.  The  elements  of  v  are  required  to  be  white  with  zero-mean,  Gaussian  statis¬ 
tics  described  by 


£{v(t,)}  =  o 

(2.22) 

~  TL(ti)  ti  =  tj 

E{v(ti)vr(tJ)}  = 

(2.23) 

0  t{  ^  tj 
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where  R  is  symmetric,  positive  definite,  and  may  be  interpreted  as  the  covariance  of  v 
at  time  t;.  Furthermore,  the  measurement  noises  for  each  measurement  type  are  usually 
assumed  uncorrelated  with  each  other.  This  assumption  results  in  a  diagonal  R  matrix. 
When  R  is  diagonal,  the  measurements  may  be  incorporated  sequentially  at  each  measure¬ 
ment  time.  For  a  standard  linear  Kalman  filter,  the  final  result  does  not  depend  on  the 
order  in  which  the  measurements  are  incorporated. 

The  filter  state  estimate  x  is  interpreted  as  the  mean  value  of  a  Gaussian  conditional 
probability  density  function  (PDF)  with  covariance  P.  The  Gaussian  requirement  implies 
that  x  is  also  the  mode  (most  likely  value)  for  the  conditional  PDF.  This  PDF  is  condi¬ 
tioned  on  ,  the  measurement  history  through  time  t{.  This  conditional  PDF  is  described 
by 

/x,e.„ Zm,(*  I  Z-)  =  [(2*)“/3  I  P(<,+  )  lV2]  "'exp  j-i  [*  -  X(t,+  )]T  P(Crl  [<  -  x(*t)]  }  (2.24) 


(Superscripts  and  appearing  with  U  are  used  throughout  this  thesis  to  indicate 
“just  after”  and  “just  before”  a  measurement  update,  respectively).  The  function  of  the 
Kalman  filter  is  to  propagate  this  conditional  probability  density  function  forward  in  time 
between  measurements  and  to  update  the  density  function  when  measurements  are  avail¬ 
able.  The  filtering  process  is  started  from  the  initial  estimates  of  the  state  vector,  Xo,  and 
covariance,  Po,  where 


x0  =  £{x(t0)}  (2.25) 

P0  =  £{[x(t0)  -  xo][x(t0)  -  xo]r}  (2.26) 

The  state  estimate  and  covariance  are  propagated  from  measurement  time  tf  (or  initial 
time)  to  measurement  time  tt“+1  by  (numerical)  integration  of 

x(t/ti)  =  F(t)x(t/ti)  +  B(t)u(t)  (2.27) 

P  (t/U)  =  F(t)P(t/U)  +  P{t/ti)FT(t)  +  G(t)Q(t)GT(t)  (2.28) 
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The  notation  “t/t,”  is  used  to  indicate  time  t  in  the  interval  |t+  where  till  mea¬ 

surements  through  time  t{  have  been  incorporated  in  the  filter  estimates  x  and  P.  At 
measurement  times  U,  these  two  quantities  are  updated  using 


K  (*,) 

=  P(t-)HT(t,)  [H(tt)P(f-)Hr(t,)  +  R(tt)]_1 

(2.29) 

x(C) 

=  x(f“ )  +  K(fj)  [z,-  -  H(f,)x(£“ )] 

(2.30) 

P(*,+  ) 

=  P(f-)-K(f,)H(t,)P(tr) 

(2.31) 

The  quantity  K  is  called  the  Kalman  filter  gain.  Note  that  K  and  P  do  not  depend  on 
the  history  of  the  measurement  vectors  z1(  Z2,  Z3 ,...Zj.  This  important  fact  allows  the  filter 
gain  and  covariance  to  be  precomputed.  However,  this  holds  true  only  for  the  case  of  linear 
dynamics  and  measurements.  The  derivation  of  these  equations  is  presented  in  Chapter  5 
of  reference  [8].  Differences  between  the  standard  Kalman  filter  structure  presented  here 
and  the  Kalman  filter  actually  implemented  for  this  research  are  discussed  later. 

2.4  Time  Correlated  Errors  and  Residual  Characteristics 

The  basic  Kalman  filter  presented  above  assumes  all  driving  and  measurement  noises 
are  white.  The  plot  of  the  power  spectral  density  (PSD)  of  such  a  white  noise  is  a  hori¬ 
zontal  line  over  all  frequencies;  a  white  noise  contains  equal  power  at  all  frequencies.  A 
continuous-time  white  noise  is  therefore  an  infinite  power  process  and  thus  cannot  really 
exist.  “Whiteness”  is  still  a  valid  and  useful  characteristic  for  many  noise  models  provided 
the  white  noise  adequately  models  the  real  noise  within  the  bandpass  of  the  system  of 
interest.  In  cases  where  the  time  correlation  of  th.  noise  cannot  be  neglected,  the  concept 
of  a  shaping  filter  is  used  to  model  the  time  correlation  properties  of  the  noise.  Such  a 
noise  is  modeled  by  one  or  more  additional  states,  xj,  which  are  augmented  to  the  basic 
state  vector.  The  Kalman  filter  described  in  this  thesis  involves  three  simple  types  of 
shaping  filters:  the  random  bias,  the  random  walk  (Brownian  motion),  and  the  first-order 
Gauss-Markov  process. 

The  term  random  bias  describes  a  quantity  with  an  unknown,  random,  initial  value 
with  a  Gaussian  PDF.  The  value  of  the  quantity  does  not  change  after  the  initial  time;  it 
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remains  a  constant  bias.  The  differential  equation  describing  this  concept  is  simply 

x(t)  =  0  (2.32) 

with 

x(t0)  =  x0,  x0  £  N[0,<r!j 

The  term  random  walk  is  a  synonym  for  Brownian  motion.  Since  white  noise  is  con¬ 
ceptually  the  derivative  of  a  Brownian  motion  process,  the  differential  equation  describing 
random  walk  is 

x(t)  =  w(f)  (2.33) 

Quantities  that  are  known  to  be  constant  biases  are  frequently  modeled,  in  the  Kalman 
filter  design,  with  random  walk  shaping  filters.  In  these  cases,  the  “w(t)”  appearing  in 
Equation  (2.33)  is  a  small  strength  pseudonoise  employed  simply  to  prevent  the  corre¬ 
sponding  elements  of  the  covariance  matrix  and,  in  turn,  the  Kalman  gain  matrix,  from 
becoming  and  remaining  zero.  If  the  filter  covariance  associated  with  a  particular  state 
becomes  zero,  then  subsequent  measurements  will  have  no  impact  on  the  filter  estimate  of 
that  state. 

The  term  first-order  Gauss-Markov  describes  a  process  with  exponentially  decaying 
time  correlation.  The  differential  equation  describing  such  a  process  is 

x(t)  =  -^x(t)  +  w(t)  (2.34) 

where  r  is  the  correlation  time.  The  autocorrelation  function  and  PSD  plots  for  a  first- 
order  Markov  process  are  shown  in  Figure  2.4.  The  first-order  Markov  process  shaping 
filter  is  seen  to  be  a  simple  low-pass  filter;  it  attenuates  the  high-frequency  content  of 
the  white  noise.  More  complex  (higher  order)  shaping  filters  may  be  synthesized  to  match 
the  PSD  of  almost  any  real  noise,  but  no  such  filters  were  considered  necessary  for  this 
research. 

The  quantity  [z<  -  Hx]  appearing  in  Equation  (2.30)  is  called  the  measurement  resid¬ 
ual,  rj.  The  residual  vector  is  the  difference  between  the  actual  measurement  vector  z,,  at 
time  f,  and  the  “expected"  measurement  vector,  H(ti)x(t“).  That  is, 

r,  s  r(tj)  =  z,  -  Hx  (2.35) 
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Ideally,  the  statistics  of  the  residual  vector  residuals  are  described  by 

E{n}  =  0  (2.36) 

E{rirf}  =  H(MP(*r)Hr(tO  +  R(ti)  (2.37) 

An  example  of  such  an  ideal  residual  time  history  is  shown  in  Figure  2.5.  When  all 
significant  processes,  including  time  correlated  noises,  are  adequately  described  in  a  state- 
space  model,  and  the  Kalman  filter  based  on  this  model  has  is  properly  tuned ,  the  true 
estimation  error  e(t)  defined  by 

e(t)  =  x(t)  —  x(t)  (2.38) 

can  be  shown  to  have  Gaussian  statistics  described  by 

£{e(t)>  =  0  (2.39) 

£{e(t)er(t)}  =  P(t)  (2.40) 

Thus  a  correctly  tuned  Kalman  filter  based  on  the  complete  state-space  model  is  an  unbi¬ 
ased  estimator;  its  estimates  have  zero  mean  error.  The  state  covariance  estimate,  P(f), 
propagated  and  updated  by  the  filter  is  seen  to  be  equal  to  the  covariance  of  the  true 
estimation  error. 

2.5  Extended  Kalman  Filter  and  Error  State  Concepts 

In  many  cases,  the  dynamics  and/or  the  measurements  involved  in  a  state-variable 
model  cannot  be  adequately  described  by  linear  equations.  In  such  cases,  the  nonlinear 
equations  corresponding  to  Equations  (2.18)  and  (2.21)  are 


x(t)  = 

f(x(t),  u(<),f]  +  G(t)w(t) 

(2.41) 

z(ti)  = 

h[x(ti),ti]  +  v(ti) 

(2.42) 

where  f  is  the  nonlinear  dynamics  function  and  h  is  the  nonlinear  measurement  function. 
Note  that  both  the  dynamics  driving  noise  and  the  measurement  noise  are  still  linearly 
additive. 
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MEASUREMENT  RESIDUAL 


Figure  2.4.  First-Order  Markov  Characteristics 
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Figure  2.5.  Example  of  a  “White”  Residual  Sequence 
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The  extended  Kalman  filter  provides  a  state  estimation  algorithm  for  such  nonlinear 
systems.  It  can  be  derived  from  the  standard  Kalman  filter  using  a  linear  perturbation 
procedure  based  on  first-order  Taylor  series  approximations  [9].  The  state  and  covariance 
propagation  from  time  £»  to  time  £t+i  is  accomplished  via  (numerical)  integration  of  the 
equations 


x(t/ti) 

=  f[x( £/*,•),  u(t),t] 

(2.43) 

p  (t/ti) 

=  F[t;  x(t/ti)]P{t/ti)  +  P(£/ti)Fr[£;x(t/ti)]  +  G(t)Q(t)Gr(t) 

(2.44) 

where 

0X  X=X(t/t<) 

(2.45) 

At  a  measurement  time  £,,  the  extended  Kalman  filter  updates  the  state  and  covariance 

using  the  equations 

K(t«)  = 

P (t~ )  Hr[ti;  x(t" )]  {Hftjj  x(t- )]  P(t- )  H r[ti;  x(tf )]  +  R(ti )}-' 

(2.46) 

x(tf)  = 

x(t~)  +  K (ti)  fa  -  h[x( £*),*«]} 

(2.47) 

P  (tf)  = 

p(tn-K(t,)H[t<;x(tn]P(er) 

(2.48) 

where 

H[£i;x(£-)j= 

OX  X=X(t~) 

(2.49) 

Note  the  structural  similarity  of  the  extended  Kalman  filter  equations  to  the  stan¬ 
dard  Kalman  filter  equations  presented  earlier.  The  major  differences  are  that  the  non¬ 
linear  dynamics  equation  is  used  for  propagation  of  the  state  estimates  and  the  nonlinear 
measurement  equation  is  used  in  formulating  the  measurement  residual.  The  matrices  F 
and  H  are,  in  general,  functions  of  the  state  estimates  as  a  result  of  the  linearization  pro¬ 
cess.  This  state  dependency  prevents  the  filter  gain  and  covariance  estimates  from  being 
precomputed.  The  accuracy  and  stability  of  an  extended  Kalman  filter’s  state  and  covari¬ 
ance  estimates  are  dependent  on  the  adequacy  of  the  linearization  process  for  the  system 
of  interest.  As  a  result  of  the  the  approximations  involved  in  the  calculation  of  F  and 
H,  the  state  estimates  produced  by  an  extended  Kalman  filter  cannot  be  mathematically 
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proven  to  be  optimal.  Even  so,  extended  Kalman  filters  may  work  very  well  for  certain 
applications  [9] . 

Another  important  concept  used  in  this  thesis  is  that  of  error  states.  Instead  of 
directly  estimating  the  numeric  value  of  some  state,  it  may  be  more  convenient  to  estimate 
the  error,  Sx,  in  the  true  value  of  that  state  relative  to  some  nominal  value  xn.  That  is, 

£x(t)  =  x(f)  -  xn(t)  (2.50) 

and  the  Kalman  filter  estimate  of  Sx  is  Sx.  Since  the  nominal  value  is  known  for  all  t  €  T, 
the  (indirect)  filter  estimate  of  the  whole  valued  state  is  calculated  as 

x{t)  =  xn(t)  +  8x{t)  (2.51) 


2.6  Simulation  Software 

All  of  the  data  processing  described  in  this  thesis  is  conducted  with  the  aid  of  two 
major  software  packages,  MSOFE  [l]  and  MATRIXx[6\-  An  overview  of  the  purpose, 
capabilities,  and  use  of  each  of  these  software  packages  is  provided  below. 

MSOFE  is  the  acronym  for  the  Multimode  Simulation  for  Optimal  Filter  Evaluation 
software  package.  MSOFE  was  jointly  developed  by  the  Avionics  Laboratory  of  Wright- 
Patterson  AFB,  Ohio,  and  Integrity  Systems  Inc.  of  Winchester,  Massachusetts.  It  is 
a  multimode  simulation  tool  for  designing  and  evaluating  integrated  systems  based  on 
optimal  (Kalman)  filtering  techniques.  The  primary  modes  of  use  are  Monte  Carlo  analysis 
and  covariance  analysis.  These  two  modes  may  be  used  separately  or  simultaneously.  The 
conceptual  and  mathematical  basis  of  these  two  analysis  techniques  is  described  in  reference 
[8:325-341].  MSOFE  is  written  in  ANSII  standard  FORTRAN  77  and  can  be  used  on  any 
reasonably  fast  computer  for  which  there  is  a  FORTRAN  77  compiler  available. 

An  executable  MSOFE  program  consists  of  am  executive  routine,  56  fixed  subrou¬ 
tines,  and  at  least  14  user-written  subroutines.  It  is  these  14  user  subroutines  that  allow 
MSOFE  to  be  customized  to  handle  almost  any  Kalman  filter  problem.  The  user  sub¬ 
routines  required  for  this  research  involved  a  significant  amount  of  original  programming. 
Although  MSOFE  was  conceived  as  a  simulation  tool  in  which  the  simulated  measure¬ 
ments  occur  at  synchronous  intervals,  it  is  easily  modified  to  filter  empirically  recorded 
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data  where  the  real  measurements  occur  asynchronously.  MSOFE  has  two  companion  pro¬ 
grams  called  PROFGEN  and  MPLOT.  PROFGEN  generates  simulated  trajectory  data  for 
use  in  MSOFE  simulations  while  MPLOT  performs  statistical  calculations  and  provides 
high  quality  plots  of  MSOFE  and  PROFGEN  output.  PROFGEN  was  not  required  in 
this  research  and  MPLOT  was  not  used  because  it  requires  commercial  plotting  routines 
which  are  not  available  on  AFIT  computers.  More  information  on  the  use  and  capabilities 
of  MSOFE  is  available  in  the  MSOFE  User’s  Manual  [1]. 

MATRIXx  is  a  commercial  software  product  developed  by  Integrated  Systems  Inc. 
of  Palo  Alto,  California  [6j.  It  is  a  powerful,  general  purpose,  interactive  program  for 
the  computer-aided  design  and  analysis  of  control  systems.  MATRIXx  implements  an 
interpreted  programming  language  similar  to  FORTRAN  and  it  can  read  data  files  created 
by  FORTRAN  programs.  Only  a  small  portion  of  its  capabilities  were  used  in  thisthesis 
effort.  MATRIXx  was  used  for  interpolation  and  plotting  of  the  MSOFE  output.  Several 
small  FORTRAN  and  MATRIXx  programs  are  used  to  handle  routine  tasks  such  as  reading 
CIRIS  data  tapes,  producing  correctly  formatted  input  data  files  for  the  MSOFE  runs,  and 
plotting  data. 

2. 7  Summary 

This  chapter  introduces  the  spatial  reference  frames,  Kalman  filter  theory,  and  sim¬ 
ulation  software  used  in  this  research.  The  primary  function  of  CIRIS  is  to  estimate 
position  and  velocity  quantities.  The  WGS  84  geodetic  reference  system  is  chosen  as  the 
basic  reference  frame  for  defining  position  and  velocity.  The  WGS-84  system  defines  the 
earth-centered  earth-fixed  (E-frame)  rectangular  coordinate  frame  and  the  relationship 
between  E-frame  coordinates  and  geodetic  latitude,  longitude,  and  altitude  coordinates. 
Several  additional  reference  frames  are  required  to  simplify  development  and  specification 
of  the  Advanced  CIRIS  filter.  These  include  the  north — west — up  local-level  navigation 
frame  ( N -frame),  the  true  frame  (T-frame),  the  platform  frame  (E-frame),  the  computer 
frame  (C-frame),  and  the  body  frame  (E-frame).  These  frames  and  the  transformations 
between  them  are  specified  in  Sections  2.1  and  2.2. 

The  concept  and  mathematical  structure  of  the  standard  linear  Kalman  filter  is 
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discussed  in  Section  2.3.  This  discussion  is  augmented  with  the  concepts  of  time  correlated 
errors,  ideal  measurement  residual  characteristics,  the  extended  Kalman  filter,  and  error 
states  in  Sections  2.4  and  2.5. 

The  results  presented  in  Chapters  IV  and  V  sire  obtained  with  the  aid  of  two  major 
software  packages.  The  first  is  the  Multimode  Simulation  for  Optimal  Filter  Evaluation 
(MSOFE)  program.  Standard  MSOFE  subroutines  implement  the  propagation  and  mea¬ 
surement  processing  procedures  generic  to  any  Kalman  filter  application  while  additional, 
user-written,  subroutines  customize  MSOFE  for  a  particular  application.  The  second  soft¬ 
ware  package  is  MATRIXx ,  a  control  system  design  and  analysis  program.  MATRIXx  is 
used  here  only  for  data  interpolation  and  data  plotting.  These  two  software  packages  are 
described  in  more  detail  in  Section  2.6 


2-17 


III.  70-State  Filter  Description 


The  Advanced  CIRIS  filter  described  in  this  chapter  is  basically  identical  to  the  70- 
state  filter  described  by  Solomon  [13]  except  for  four  significant  changes.  The  first  change 
is  that  the  measurement  models  are  developed  using  an  extended  Kalman  filter  formulation 
rather  than  a  linearized  Kalman  filter  formulation.  The  second  change  is  that  the  filter  is 
updated  with  barometric  altitude  measurements  in  addition  to  the  transponder  range  and 
range-rate  measurements.  The  third  change  is  that  the  atmospheric  propagation  delay 
is  modeled  as  being  range-dependent  rather  than  range-independent,  and  that  the  atmo¬ 
spheric  delays  associated  with  the  different  transponders  have  non-zero  correlation.  The 
fourth  change  is  the  compensation  for  the  position  and  motion  of  the  CIRIS  interrogator 
antenna  relative  to  the  INS.  In  the  past,  this  lever  arm  correction  was  provided  by  the 
CIRIS  II  software. 

3. 1  Notation 

The  primary  function  of  the  filter  is  to  produce  accurate  estimates  of  the  position 
and  velocity  of  the  CIRIS  LN-39  INS  during  a  specified  time  interval.  The  mathematical 
equations  which  describe  the  filter’s  computations  require  rather  complex  notation.  This 
section  defines  terms  and  explains  notation. 

The  measurements  processed  by  the  filter  are  range  and  range-rate  measurements 
provided  by  the  Cubic  Range  and  Range-Rate  Subsystem  (RRS)  and  a  barometric  altitude 
measurement  provided  by  the  central  air  data  computer  (CADC).  The  principle  compo¬ 
nents  of  the  Cubic  RRS  are  the  transponder  electronics,  transponder  antennas,  interrogator 
electronics,  and  interrogator  antenna.  The  filter  estimation  of  the  position  and  velocity  of 
the  INS  requires  estimates  of  the  position  and  velocity  of  the  transponder  antennas  and 
the  interrogator  antenna. 

Let  the  position  vectors  of  the  jth  transponder  antenna,  the  interrogator  antenna, 
and  INS  be  denoted  by  Pt,,  Pa,  and  Pj,  respectively.  In  the  .E-frame,  these  positions 
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are  shown  in  terms  of  their  components  as 

kt  r*r  f  x,  r 


The  familiar  “dot”  (')  notation  is  used  to  indicate  differentiation  with  respect  to  time. 
Thus  the  .E-frame  velocities  of  these  three  points  are  denoted  as 


The  “tilde”  (“)  symbol  is  used  to  indicate  measured  or  calculated  quantities  where  the  mea¬ 
surement  or  calculation  process  is  assumed  to  introduce  error.  For  example,  transponder 
positions  are  provided  to  the  filter  from  a  database  containing  the  surveyed  (measured) 
transponder  WGS  84  lathd  h,  longitude,  and  altitude.  The  INS  calculates  its  own  position 
based  on  its  initializati  ,u  and  subsequent  inertial  measurements.  The  interrogator  antenna 
position  is  calculated  relative  to  the  INS  position  assuming  a  fixed  5-frame  offset  vector 
and  knowledge  of  the  aircraft  orientation.  These  measured  or  calculated  quantities  are 
denoted  as 


Error  (difference)  quantities  sire  indicated  with  a  delta  (<5)  prefix.  For  example,  the  error 
in  the  indicated  INS  position  relative  to  the  true  INS  position  is  denoted  as 

q  E  p  ^  E  *  E 

Xj  X/ 

SY i  =7/  -  Vi 

SZj  Zj  Zj 

Finally,  filter  estimates  of  a  quantity  are  indicated  by  a  “hat”  (~)  symbol.  For  example,  if 
the  E-frame  X ,  7,  and  Z  errors  in  the  INS  position  are  three  of  the  filter  states,  the  filter 
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estimates  of  these  states  are  denoted  as 


6X1 
SYi 
6Z / 

The  (indirect)  filter  estimate  of  the  true  INS  position  is  calculated  by  subtracting  the  filter 
estimate  of  the  error  in  the  INS  indicated  position  from  the  INS  indicated  position.  This 
concept  is  formulated  as 


’  Xj 

E 

'  Xj 

E 

8X1 

?I 

= 

Yi 

- 

6Y1 

Zi 

.  ^ 

SZi 

The  symbolic  order  of  operations  is  important.  For  example,  8Zi  indicates  the  filter 
estimate  of  the  error  in  the  indicated  INS  velocity  along  the  Z  axis  while  8Zj  indicates  the 
time  rate  of  change  of  the  filter  estimate  of  the  error  in  indicated  INS  position  along  the 
Z  axis. 

3.2  Error  State  Vector 

This  version  of  the  Advanced  CIRIS  filter  models  70  error  quantities.  The  error 
state  vector  is  partitioned  into  four  subvectors  which  group  related  error  states.  The 
total  error  state  vector  is  denoted  as  6x.  The  subvectors  are  defined  as  follows:  5xj 
contains  the  13  general  INS  dynamics  errors  such  as  position,  attitude,  velocity,  and  vertical 
channel  errors;  €x 2  contains  12  additional  INS  gyro  and  accelerometer  errors;  6x3  contains 
3  baro-altimeter  errors;  and  5x4  contains  all  42  transponder  related  errors.  While  many 
of  the  individual  elements  of  6x  are  defined  along  a  specific  axis  of  either  the  T-frame 
or  the  i?-frame,  others  are  defined  along  time  varying  line-of-sight  directions  or  do  not 
correspond  to  a  spatial  direction  at  all.  Thus  it  is  not  correct  to  think  of  the  whole  6x 
vector  as  being  coordinated  in  a  particular  3-dimensional  spatial  reference  frame;  the 
elements  of  6x  are  the  coordinates  of  a  point  in  a  70-dimensional  error  state-space.  In 
general,  the  elements  of  subvectors  5xi,  Sx2,  and  6x3  are  defined  relative  to  the  T-frame 
while  the  elements  of  subvector  6x4  are  defined  relative  to  the  £-frame.  The  individual 
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elements  of  6x  and  their  associated  initial  ler  value  are  defined  in  Appendix  A,  Tables 
A.l  thru  A. 6.  The  most  significant  error  states,  those  in  subvectors  6xi,6x3,  and  6x4,  are 
described  in  more  detail  below.  The  elements  of  error  state  subvector  6x2,  while  necessary 
for  the  completeness  of  the  INS  error  model,  are  more  difficult  to  visualize.  For  more 
information  on  these  error  quantities  the  reader  is  referred  to  the  Litton  LN-39  Systems 
Engineering  Analysis  Report  [14]. 

3.2.1  INS  Primary  Error  State  Subvector  The  errors  in  the  INS  indicated  position 
and  platform  wander  angle  are  modeled  by  the  four  error  states  S&y,  8Qz,  and 

6h.  Error  states  6Qx >  and  SQz  are  the  error  angles  appearing  in  Equation  (2.6). 

Together,  they  define  the  INS  latitude,  longitude,  and  platform  wander  angle  errors,  while 
error  state  Sh  defines  the  INS  altitude  (vertical  position)  error.  Their  relation  of  the  60  er¬ 
ror  angles  to  the  more  familiar  latitude  ( SL ),  longitude  (5A),  and  wander  angle  (£a)  errors  is 


and 


6L 

sin  a 

cos  a  0 

SQx  ' 

6\ 

= 

cos  a  sec  L 

-  sin  a  sec  L  0 

se  y 

6h 

0 
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6a  =  -(cos  a  tan  L)  6Qx  +  (sin  a  tan  L)6@y  +  8&z  ( 3.2 ) 


These  equations  are  taken  from  the  Litton  LN-39  Systems  Engineering  Analysis  Report 
[14].  The  error  quantity  6a  is  explicitly  computed  to  correct  the  INS  indicated  wander  angle 
for  use  in  calculating  the  state  dynamics  matrix  and  various  coordinate  transformation 
matrices. 


The  latitude  and  longitude  errors  shown  in  Equation  (3.1)  are  in  angular  units  (radi¬ 
ans).  The  transformation  which  converts  these  errors  to  linear  distance  units  coordinatized 
in  the  N-frame  is 
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(3.3) 


3-4 


where 

R  =  A(1  - /sin2  I)  (3.4) 

is  the  radius  from  the  center  of  the  WGS  84  ellipsoid  to  a  point  on  its  surface  at  latitude 
L.  When  Equations  (3.1)  and  (3.3)  are  combined,  the  result  is 


8N 
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(IZ-|-/i)sinQ  (R  +  h)  cos  a  0 
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0  0  1 

8h 

The  transformation  matrix,  Tjf,  appearing  in  Equation  (3.5)  is  used  in  the  similarity 
transformation  for  transforming  the  filter  computed  covariances  associated  with  8&x ,  8Qy , 
and  8h  into  the  corresponding  N -frame  covariances. 

The  error  in  the  INS  platform  orientation  relative  to  the  true  local-level  plane  is 
modeled  by  the  error  states  <f>x,  <t>Y ,  and  <f>Z-  These  three  error  states  are  the  misalignment 
angles  appearing  in  Equation  (2.5).  Together,  they  describe  the  platform  tilt  relative  to 
the  T-frame. 


The  error  in  the  INS  indicated  velocity  is  modeled  by  the  error  states  8Vx ,  8Vy,  and 
8Vz  .  These  three  error  states  are  defined  as  the  components  of  the  difference  vector  between 
the  INS  computed  (C-frame)  velocity  and  the  true  (T-frame)  velocity.  These  three  error 
states  may  be  transformed  into  equivalent  N -frame  velocity  errors  by  compensating  for 
the  platform  wander  angle.  The  transformation  is 


(3.6) 
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The  transformation  matrix,  Cp,  appearing  in  Equation  (3.6)  is  used  in  the  similarity 
transformation  for  transforming  the  filter  computed  covariances  associated  with  8Vx ,  8Vy , 
and  8Vz  into  the  corresponding  N -frame  covariances. 
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The  three  error  states  Shi,  SS 3,  and  8S 4  model  errors  in  quantities  used  internally  by 
the  INS  for  calculating  vertical  channel  aiding  gains.  The  LN-39  requires  baro-altimeter 
aiding  for  stabilization  of  the  vertical  channel.  Figure  3.1  is  taken  from  reference  [14]  to 
illustrate  the  LN-39  vertical  channel  error  model.  The  INS  indicated  altitude  is  lagged 
(held)  for  one  second  before  being  compared  to  the  baro-altimeter  indicated  altitude.  This 
one  second  wait  is  used  to  compensate  for  the  delay  inherent  in  baro-altimeter  operation. 
The  difference  signal,  h(t,_i )  -  hb(U),  is  used  by  four  different  vertical  channel  aiding  loops 
within  the  LN-39.  The  aiding  gains  for  these  loops  are  designated  K\,  K3,  K3,  and  A'4. 
They  are  nonlinear  functions  of  the  rate  of  change  of  the  baro-altimeter  indicated  altitude, 
hty.  The  Litton  LN-39  Systems  Engineering  Analysis  Report  defines  these  gains  as 


where 
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'  A0  -  8  if  A0  >  A  and  A0  >  38 

30  ft/ sec  otherwise  (and  initially) 


100 


1  + 


21 


(3.11) 


(3.12) 


(3.13) 


Note  that  Ao  changes  in  discrete  increments  (0,+8,  or  -8)  each  time  it  is  recalculated. 
This  implies  the  aiding  gains  also  change  in  discrete  increments. 


The  dynamics  matrix  elements  for  error  states  Shi,  SS 3,  and  6S4  are  listed  in  Ap¬ 
pendix  A,  Tables  A. 7  through  A. 9.  Note  that  the  differential  equation  for  SS 4  contains 
discontinuities  due  to  its  dependence  on  A'4.  Within  the  filter  software,  the  aiding  gains  are 
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Figure  3.1.  LN-39  INS  Vertical  Channel  Error  Model 
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recalculated  only  after  each  measurement  update  cycle  so  that  the  SS4  differential  equation 
is  continuous  throughout  each  propagation  cycle. 

The  initial  l<r  values  for  the  13  states  in  subvector  £xi  are  specified  in  Appendix  A, 
Table  A.l. 

3.2.2  Baro-Altimeter  Error  State  Subvector  The  baro-altimeter  error  is  modeled 
by  the  three  error  states  Sh\, lt  Shi, 2,  and  Sh^.  State  6h\,\  represents  the  time-correlated 
component  of  the  baro-altimeter  error.  It  is  modeled  with  a  first-order  Markov  shaping 
filter.  State  Sh^  represents  the  constant  bias  in  the  baro-altimeter  output.  It  is  modeled 
with  a  random  bias  shaping  filter.  State  Sh\,z  represents  the  altitude  scale  factor  for  the 
altitude  dependent  bias  in  the  baro-altimeter  output.  It  is  also  modeled  with  a  random 
bias  shaping  filter.  The  total  baro-altimeter  error,  Shi,,  is  defined  in  terms  of  these  three 
error  states  and  the  true  altitude  h  as 

Shb(t)  =  hf,(t)  —  h(t)  =  Shf,i(t)  -+■  Shi,2  +  h(t)  Sh (,3  (3.14) 

The  initial  l<r  values  for  these  three  states  are  specified  in  Appendix  A,  Tables  A. 3.  The 
dynamics  matrix  elements  are  listed  in  Appendix  A,  Tables  A. 7  through  A. 9. 

3.2.3  Transponder  Error  State  Subvector  The  transponder  related  errors  are  di¬ 
vided  into  those  common  to  all  transponders  and  those  concerning  individual  transponders. 
The  common  element  in  the  range  and  range-rate  measurements  from  all  transponders  is 
the  interrogator.  The  error  common  to  all  range  measurements  due  to  an  interrogator 
calibration  bias  is  represented  by  error  state  SErc •  The  error  common  to  all  range-rate 
measurements  due  to  an  interrogator  calibration  bias  is  represented  by  error  state  SErC. 
Both  error  states  are  modeled  with  random  bias  shaping  filters. 

The  errors  associated  with  an  individual  transponder  include  survey  (position)  errors 
and  a  range  measurement  error  caused  by  atmospheric  propagation  delay  along  the  line- 
of-sight  between  the  transponder  and  the  interrogator  antennas.  The  survey  errors  for 
transponder  Tj  are  represented  by  error  states  SXj, ,  SYpj ,  and  SZp,  which  model  the 
components  of  the  position  error  along  the  Xe,  Ye,  and  Ze  axes,  respectively.  These  three 
error  states  are  each  modeled  with  random  bias  shaping  filters.  The  range  error  due  to 
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the  atmospheric  propagation  delay  is  represented  by  error  state  6  At, .  This  error  state  is 
interpreted  as  a  scale  factor;  the  range  error  is  the  product  of  the  scale  factor  and  the  range. 
Error  state  6  At,  is  modeled  with  a  first-order  Markov  process  shaping  filter.  These  four 
types  of  error  are  illustrated  in  Figure  3.2.  The  driving  noises  for  the  6 At,  error  states  are 
modeled  as  being  correlated  with  each  other  with  a  positive  correlation  coefficient.  This 
additional  tuning  parameter  provides  a  crude  model  for  the  commonality  of  atmospheric 
conditions  along  the  lines-of-sight  to  each  transponder.  The  initial  Icr  values  for  these  42 
error  states  are  specified  in  Appendix  A,  Tables  A. 4  through  A. 6.  The  dynamics  matrix 
elements  are  listed  in  Appendix  A,  Tables  A. 7  through  A. 9. 

The  filter  model  includes  a  SXt,,  SYt,,  6Zt,,  and  6 At,  error  state  group  for  each  of 
ten  different  transponders  (j  =  1,2, .. .,  10).  However,  a  typical  CIRIS  test  flight  involves 
more  than  ten  different  transponders.  The  filter  software  includes  a  switching  subrou¬ 
tine  which  “timeshares”  each  of  the  transponder  error  state  groups  between  two  different 
transponders.  This  permits  the  filter  software  to  process  measurements  from  up  to  twenty 
different  transponders.  At  any  given  time,  one  member  of  each  transponder  pair  is  “active” 
(owns  the  error  state  group)  while  the  other  is  “inactive.”  If  the  filter  receives  a  measure¬ 
ment  from  the  inactive  transponder,  the  switching  subroutine  saves  the  four  error  state 
values  associated  with  the  active  transponder  and  reloads  the  four  error  state  values  previ¬ 
ously  stored  for  the  inactive  transponder.  It  then  clears  the  four  rows  and  four  columns  of 
the  covariance  matrix  associated  with  these  error  states,  saving  only  the  diagonal  elements, 
and  resets  the  four  diagonal  elements  in  these  rows/columns  to  their  previously  stored  val¬ 
ues.  The  inactive  transponder  is  thus  activated  while  the  previously  active  transponder  is 
deactivated.  This  procedure  doubles  the  number  of  allowable  transponders  at  the  expense 
of  losing  covariance  information.  The  off-diagonal  covariance  information  cannot  be  stored 
and  later  reloaded  because  any  of  the  ten  transponders  pairs  may  “switch”  at  any  time, 
thus  rendering  the  cross-correlation  information  invalid.  The  frequency  of  switching  is 
minimized  by  pairing  widely  separated  transponders  so  that  it  is  unlikely  that  the  inter¬ 
rogator  will  request  measurements  form  both  members  in  any  short  (10-30  minutes)  time 
interval. 
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Figure  3.2.  Relative  Positions  of  Transponder,  Antenna,  and  INS 
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3.3  Error  State  Stochastic  Differential  Equation 

The  dynamics  of  the  error  state  vector  are  described  by  the  first-order,  linear,  piece- 
wise  continuous,  matrix  differential  equation 

6x(t)  =  F[t;  x(t)]  6x(t)  +  w(t)  (3.15) 

This  equation  is  piecewise  continuous  because  the  elements  of  the  F  matrix  that  depend 
on  the  baro-altimeter  aiding  gains  experience  discontinuities  at  measurement  times.  The 
elements  of  the  vector  w  represent  the  zero-mean  white  noises,  if  any,  driving  each  state. 
The  strengths  of  these  noises  are  represented  by  a  matrix  Q  which  is  defined  in  accordance 
with  Equation  (2.20). 

It  is  useful  to  rewrite  Equation  (3.15)  to  illustrate  the  interdependencies  of  the  error 
state  vector  partitions.  The  partitioned  form  of  Equation  (3.15)  is 

Sx  x  Fj.i  Fii2  Fii3  0  Sxi  wi 

6x2  0  F2  2  0  0  6x2  w2 

=  +  (3.16) 

6x3  0  0  F3,3  0  6x 3  w3 

6x 4  0  0  0  F4i4  6x4  w4 

This  form  of  the  dynamics  matrix  shows  that  subvectors  6x2,  6x3,  and  6x4  do  not  interact 
during  the  filter  propagation  cycle.  Cross-coupling  of  these  subvectors  occurs  only  during 
the  measurement  update  procedure.  The  non-zero  elements  of  submatrices  Fi,i,  Fj  2, 
Fi,3,  F2]2,  and  F3|3  are  given  in  Appendix  A,  Tables  A. 7  through  A. 9.  The  elements  of  the 
submatrices  Q3,  Q3,  and  Q4,  which  correspond  to  the  four  partitions  of  the  w  vector, 
are  given  in  Appendix  A,  Table  A. 10.  Pseudonoise  is  not  used  with  any  of  the  random  bias 
shaping  filter  states  because  analysis  using  real  data  indicated  the  covariances  associated 
with  these  states  remained  non-zero  throughout  the  required  time-interval.  In  general, 
pseudonoise  with  an  appropriately  chosen  strength  should  be  used  with  random  bias  errors. 

3-4  Measurement  Models 

The  filter  state  estimates  are  periodically  updated  with  three  types  of  measurements. 
A  transponder  interrogation  occurs  approximately  once  per  second;  a  successful  interroga¬ 
tion  provides  a  range  and  a  range-rate  measurement  from  one  of  the  transponders  in  the 
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current  transponder  set.  A  barometric  altitude  measurement  is  made  at  the  same  time 
as  each  transponder  interrogation.  Therefore  the  measurement  vector  z  can  be  written  in 
terms  of  its  three  components  as 

Zr(ti) 

Zr(ti)  .  (3.17) 

2hk  ( tj  ) 

where  the  three  components  are,  from  top  to  bottom,  the  transponder  range  measure¬ 
ment,  the  transponder  range-rate  measurement,  and  the  barometric  altitude  measurement. 
These  three  measurements  are  processed  one  at  a  time  (scalar  update)  by  the  MSOFE  soft¬ 
ware.  The  U-D  factorization  scalar  update  procedure  used  by  MSOFE  [lj  requires  a  scalar 
equation  of  the  form  of  Equation  (2.42)  be  specified  for  each  of  the  three  measurement 
types.  Thus 

h„[x(t;),*»] 

h[x(*<),*i]=  h  4x(tiMi]  (3.18) 

hkjx(tj),  £«] 

where  hr,  h^,  and  h/,4  are  each  scalar  valued  equations.  Each  of  these  three  functions  has  a 
corresponding  l-by-70  observation  matrix  row  defined  in  accordance  with  Equation  (2.49). 

All  range  and  range-rate  measurements  are  defined  in  terms  of  E-frame  coordi¬ 
nates.  They  must  be  projected  into  the  error  state  vector  space  for  correct  updating  of  the 
INS  error  and  baro-altimeter  error  states.  This  projection  operation  transforms  the  INS 
and  baro-altimeter  elements  of  the  E- frame  observation  matrix  into  the  equivalent  error 
state-vector  space  observation  matrix.  The  form  of  this  transformation  is  presented  in 
Appendix  A,  Table  A.  11  for  the  range  measurement  and  in  Table  A.  12  for  the  range-rate 
measurement. 

3.4-1  Range  Measurement  The  range  measurement  is  provided  by  the  Cubic  RRS. 
Ideally,  the  range  measurement  is  the  line-of-sight  distance  between  the  transponder  and 
interrogator  antennas  plus  a  known  bias  caused  by  propagation  delay  in  the  cables  connect¬ 
ing  the  antennas  to  their  associated  signal  processing  electronics.  The  range  measurements 
used  in  this  thesis  research  contain  an  open-loop  correction  for  the  cable  propagation  delay 
as  provided  by  the  CIRIS  II  software  [11], 
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The  following  discussion  develops  the  range  measurement  model  used  in  the  Ad¬ 
vanced  CIRIS  Kalman  filter.  As  before,  let  the  5-frame  position  vectors  of  the  jth 
transponder  antenna,  the  interrogator  antenna,  and  INS  be  denoted  by  Pf:  ,  PA,  and 
Pf ,  respectively.  The  line-of-sight  range,  Rj,  between  the  two  antennas  is  a  nonlinear 
function  of  these  position  coordinates.  The  functional  form  of  this  Euclidian  distance  is 


R:  = 


=  y/{xT)  -  XA )2  +  (yt,  -ya)2  +  (zt,  -  zAy- 


(3.19) 


This  nonlinear  dependence  dictates  the  use  of  the  extended  Kalman  filter  measurement 
update  procedure,  Equations  (2.46)  through  (2.49). 

The  range  measurement  function,  hr,  is  defined  as  the  sum  of  the  true  line-of- 
sight  range,  Rj,  the  range  dependent  atmospheric  propagation  delay,  Rj6AT} ,  and  the 
interrogator  range  measurement  calibration  bias,  6Ere.  This  measurement  is  corrupted  by 
a  random  measurement  noise,  vr.  In  accordance  with  Equation  (2.42),  the  assumed  model 
of  the  range  measurement  provided  by  transponder  Tj  is 


zr(ti)  =  hr[tj,  x(t{)j  +  ur(tj) 


(3.20) 


where 

hr[ti,x(ti)]  =  +  5j(x(tj)]tfi4r,(<t)  +  tiEre(ti)  (3.21) 

The  range  measurement  noise  is  modeled  as  a  zero-mean  white  noise.  The  range  measure¬ 
ment  noise  strength  is  specified  in  Appendix  A,  Table  A. 14. 

The  corresponding  range  observation  vector,  Hr,  is  calculated  using  Equation  (2.49). 
The  result,  shown  here  as  the  bracketed  terms  in  the  product  Hr6x,  is 


Hr[*<;x(<j  )]$x(tj)  =  + 

Xjv  -  XA 

[ 
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6  Zt 

[«,]  6At,  +  (!)  6Er 


(3.22) 


Note  that  hr  and  Hr  both  depend  on  the  E- frame  position  of  the  interrogator  antenna, 
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but  the  position  error  states  estimated  by  the  filter  are  for  the  INS  position.  The  antenna 
location  must  be  calculated  using  the  filter  estimates  of  the  INS  position  error  and  knowl¬ 
edge  of  the  offset  (lever  arm)  between  the  INS  platform  and  the  interrogator  antenna.  The 
offset  vector  is  constant  in  the  5-frame  but  depends  on  the  aircraft  position  and  orien¬ 
tation  in  the  5-frame.  This  offset  vector,  M ,  is  illustrated  in  Figure  3.2.  It  is  defined 
by 
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(3.24) 


The  5-frame  position  coordinates  of  the  interrogator  antenna  are  estimated  by  the  filter 
software  as 


(3.25) 
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The  5-frame  components  of  M  are  determined  from  Equation  (3.24)  using  filter  estimates 
of  the  six  arguments. 


3.4.2  Range-Rate  Measurement  The  range-rate  measurement  is  also  provided  by 
the  Cubic  RRS.  Ideally,  the  range-rate  is  simply  the  magnitude  of  the  relative  velocity 
between  the  transponder  and  interrogator  antennas.  This  (5-frame)  range-rate  is  found 
by  differentiating  Equation  (3.19).  The  result  is 


R,[x(t),t]  = 


dR,  dRj  dXTj  OR,  dXA 

dt  ~  dXTi  dt  +  dXA  dt  + 


dXTj 

- 


dR,  dZf  OR,  dZA 

j - - - L  - 1 - _ 

dZr,  dt  dZA  dt 


Zt,  -  ZA 

R, 


ZA  (3.26) 
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The  range-rate  measurement  function,  h,.,  is  defined  as  the  sum  of  the  true  line- 
of-sight  range-rate,  Rj,  and  the  interrogator  range-rate  measurement  calibration  bias, 
6 Ere ■  This  measurement  is  corrupted  by  a  random  measurement  noise,  tv.  In  accordance 
with  Equation  (2.42),  the  assumed  model  of  the  range-rate  measurement  provided  by 
transponder  Tj  is 

Z  r{ti)  =  hr[ti,x{ti))  +  Vr(tj)  (3.27) 


where 


hf[ti,x(ti)]  =  flj[x(*i)]  +  SErc{ti) 


(3.28) 


The  range-rate  measurement  noise  is  modeled  as  a  zero-mean  white  noise.  The  range-rate 
measurement  noise  strength  is  specified  in  Appendix  A,  Table  A. 14. 

The  corresponding  range-rate  observation  vector,  Hr,  is  calculated  using  Equa¬ 
tion  (2.49).  The  result,  shown  here  as  the  bracketed  terms  in  the  product  Hr<5x,  is 


Hr[ti;x(t{  )]6x(*i)  = 

'(At,  -  XAfXA  +  (XT.  -  Aa)(?t,  -  ?a)Ya  +  (XT.  -  XA)(ZTi  -  ZA)ZA  -  R\  Xa 


^t, 


6X 


a/ 


(Yt,  -  Ya)(Xt,  -  XA)XA  +  (Yt.  -  YA)2YA  +  (YTj  -  Ya)(ZTj  -  ZA)ZA  -  R^  Y A 

R%. 

(ZTj  -  ZA)(XTi  -  XA)XA  +  (ZTj  -  ZA){YTl  -  Ya)Ya  +  ( Z t,  -  ZA)2ZA  -  r*za 

Rj,. 
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SZi 
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6X  i  — 

'yt,-Ya 
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(3.29) 


Note  that  hf  and  Hf  depend  on  the  E-frame  velocity  of  the  interrogator  antenna  as  well 
as  well  as  its  E-frame  position,  but  the  velocity  error  states  estimated  by  the  filter  are  for 
the  INS  velocity.  The  E-frame  velocity  of  the  interrogator  antenna  is  found  as  follows. 
Since 

[Pa}E  =  [Pi)E  +  E  Mf  (3.30) 


then 

[Pa]E  =  [Pif  +  [Mf  (3-31) 
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Also,  since 

[M}e  =  C%[M]B  (3.32) 

then 

[M)e  =  CEb[M]b  (3.33) 

because  the  5-frame  representation  of  M  is  constant.  The  derivative  of  C §  is  found  using 
the  identity 

&B  ~  ^B^EB  (3.34) 

where  f2BB  is  the  skew-symmetric  form  o{ujbb,  the  5-frame  to  5-frame  angular  velocity 
vector  expressed  in  the  5-frame.  Because  a >BB  is  not  directly  available,  it  is  calculated  as 

B  B  ,  B 

usEB  ~  ^en+Uhb 

~  Cn^en+^nb  (3.35) 


where  u is  computed  from  the  latitude  and  longitude  rates  and  wBB  is  computed  from 
the  heading,  pitch,  and  roll  rates.  These  two  angular  rotation  rate  vectors  are 


*EN  = 


A  cos  L 


A  sin  L 


4>+  4<  sin  9 
B 

<*>NB  =  9  cos  <{>  -  t/>  sin  <p  sin  9 

9  sin  <j>  +  4 >  cos  <fi  cos  9 


When  these  equations  are  combined,  and  filter  estimates  used  for  all  needed  matrix  ele¬ 
ments,  the  filter  estimate  of  the  5-frame  velocity  of  the  interrogator  antenna  is  computed 


'■E-~~B 


[PaT=  ya  =  y i  +  cBnEB  YM 


(3.36) 


3-4-3  Barometric  Altitude  Measurement  The  altitude  measurement  is  provided  by 
the  barometric  altimeter.  Ideally,  this  measurement  indicates  the  true  altitude  of  the  INS. 
In  order  to  account  for  limitations  inherent  in  the  estimation  of  altitude  based  on  baro¬ 
metric  pressure  and  other  error  sources,  the  baro-altimeter  output  is  modeled  as  the  sum 
of  the  true  altitude,  a  constant  bias  error,  a  time  correlated  bias  error,  a  scale  factor  error, 
and  a  random  measurement  noise,  t >hh.  In  order  to  eliminate  the  true  altitude  from  the 
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measurement  and  isolate  the  altitude  errors,  the  measurement  presented  to  the  filter  is  the 
difference  between  the  barometric  altitude  and  the  INS  indicated  altitude.  In  accordance 
with  Equation  (2.42),  the  assumed  model  of  the  barometric  altitude  measurement  is 


ZfcJ* «)  =  hht[t<,  x(f<)]  +  vhfi(ti)  (3.37) 

where 

=  hb(ti)-h(ti)  (3.38) 

=  [h{ti)  +  6hb{ti)]-[h(ti)  + 6h{t{)}  (3.39) 

=  Shb(tt)  -  6h(ti)  (3.40) 

=  [^/i;,i(ti)  +  Shb2  +  h(ti)6hi>$]  —  Sh  (3.41) 


The  corresponding  barometric  altitude  observation  vector,  H/,,,,  is  calculated  using  Equa¬ 
tion  (2.49).  The  result,  shown  here  as  the  bracketed  terms  in  the  product  Hh,6x,  is 

H*Jtt;x(t-)]M*<)  =  W^i»i(f«)  +  [1]^62(*,-)  +  [h(ti)\  Shb3(tt)  -  [l]6h  (3.42) 

The  baro-altimeter  output  is  sampled  within  50  milliseconds  of  each  transponder 
interrogation.  Since  it  is  not  separately  time-tagged,  the  baro-altimeter  measurement  is 
assumed  to  have  the  same  time-tag  as  the  associated  transponder  range  measurement  [11]. 

3.5  Summary 

This  chapter  presents  the  structure  of  the  70-state,  MSOFE-based,  extended  Kalman 
filter  that  is  the  current  incarnation  of  the  Advanced  CIRIS  filter.  The  (error)  state  vector 
is  partitioned  into  four  error-state  groups.  The  first  group  contains  13  INS  position, 
velocity,  platform  tilt,  and  vertical  channel  aiding  errors;  the  second  group  contains  12 
gyro  and  accelerometer  errors;  the  third  group  contains  three  baro-altimeter  errors;  and  the 
fourth  group  contains  42  transponder  calibration,  position,  find  atmospheric  propagation 
errors.  Of  the  42  transponder  states,  40  represent  ten  sets  of  the  four  error  states  associated 
with  each  of  up  to  ten  different  transponders.  The  filter  software  permits  the  use  of  up  to 
twenty  different  transponders  through  a  state  time-sharing  procedure. 
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The  filter  uses  three  types  of  measurements.  Range  and  range-rate  measurements  are 
provided  by  the  transponders,  while  barometric  altitude  is  provided  by  a  baro-altimeter. 
The  range  and  range-rate  measurement  functions  are  nonlinear.  Because  of  the  position 
offset  between  the  CIRIS  antenna  and  the  INS,  a  procedure  for  compensating  for  “lever- 
arm”  effects  is  developed. 
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IV.  70-State  Filter  Evaluation:  Ground  Test  Data 

This  chapter  describes  the  use  of  the  70-state  filter  to  process  real  measurements 
collected  during  a  slow  speed  ground  test. 

4-1  Data  Collection  Procedure 

The  data  used  in  this  analysis  was  collected  on  21  July  1989,  during  a  series  of  tests  at 
the  65 85t/l  Test  Group’s  Rocket  Sled  Test  Track.  The  CIRIS  LN-39  INS,  transponder  inter¬ 
rogator,  baro-altimeter,  and  HP-1000  computer  were  placed  in  a  Recovery  Support  Vehicle 
(RSV).  Starting  from  a  known  point  on  the  track,  the  RSV  was  driven  either  north  or  south 
along  the  track  while  the  CIRIS  equipment  recorded  the  INS  data,  barometric-altitude, 
and  the  transponder  measurements.  At  the  same  time,  the  position  of  the  RSV  was 
recorded  by  an  independent  track  data  acquisition  system  (TDAS).  Both  sets  of  recorded 
data  were  referenced  to  a  common  IRIG  time  base.  The  starting  configuration  is  illustrated 
in  Figure  4.1.  The  offset  dimensions  between  the  INS  and  the  TDAS  sensor  are  listed  in 
Table  4.1. 

This  series  of  tests  involved  six  different  transponders.  The  transponder  ID’s  and 
their  corresponding  WGS  84  coordinates  are  listed  in  Table  4.2.  The  position  of  these 
transponders  relative  to  the  track  is  shown  in  Figure  4.2.  The  initial  INS  alignment  lo¬ 
cation  and  starting  point  for  each  of  the  northward  runs  was  32°  53'  6.5648"  N  latitude, 
106°  8'  59.6288"  W  longitude,  and  3991  feet  altitude. 

The  measurements  from  transponder  211  are  not  used  because  the  range  measure¬ 
ments  contain  an  anomalous  oscillation  with  a  period  of  about  90  seconds.  The  oscillation 
is  not  apparent  in  the  associated  range-rate  measurements.  The  oscillation  is  shown  in 
Figure  4.3.  This  plot  is  obtained  by  using  the  TDAS  indicated  positions  to  calculate  the 
expected  range  to  the  transponder  at  each  point  along  the  track.  The  expected  range 
is  then  subtracted  from  the  measured  range  at  each  measurement  time.  The  bias  shown 
in  Figure  4.3  is  most  likely  due  the  combination  of  an  error  in  the  surveyed  position  of 
transponder  211  and  the  offset  between  the  CIRIS  antenna  and  the  TDAS  sensor.  However 
the  oscillation  cannot  be  explained  at  this  time.  The  oscillation  is  not  present  when  the 
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Figure  4.1.  Equipment  Configuration  for  the  CIRIS/RSV  Tests 


Table  4.1.  IKS  Position  Offsets 


South-to-North  Run 

North-to-South  Run 

23.02  ft  North 
-5.69  ft  East 

4.75  ft  Above 

-23.02  ft  North 

2.60  ft  East 

4.75  ft  Above 
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Table  4.2.  Transponder  Locations  for  the  CIRIS/RSV  Tests 


Transponder  ID 

Latitude 

Longitude 

Altitude 

005 

33°  01'  36.1472" 

-106°  08' 20.7404" 

4339  ft 

102 

32°  55'  58.5986" 

-106°  08'  50.3339" 

4074  ft 

181 

33°  44'  58.035" 

-106°  22'  14.630" 

7932  ft 

211 

33°  17'  55.999" 

-106°  31'  44.311" 

8842  ft 

212 

32°  47'  16.418" 

-105°  49' 15.474" 

9202  ft 

216 

32°  42'  12.235" 

-106°  07'  38.907" 

4481  ft 

RSV  is  stationary.  In  two  other  track  runs  in  which  the  RSV  speed  was  purposely  varied, 
the  frequency  of  the  oscillation  increased  in  proportion  to  the  speed.  In  computer  rims 
where  the  transponder  211  range  measurements  are  processed  by  the  filter,  the  oscillation 
appears  in  plots  of  the  measurement  residuals  and  in  both  the  transponder  and  INS  po¬ 
sition  error  states.  An  example  of  the  oscillating  residuals  is  shown  in  Figure  4.4.  The 
dashed  lines  in  this  plot  awe  the  filter  calculated  residual  ±1<7  bounds.  The  oscillation  is 
also  apparent  in  the  transponder  211  range  measurement  residuals  of  the  current  CIRIS  II 
filter. 


4-2  70-State  Filter  Performance 

The  original  objective  for  the  use  of  this  ground  test  data  was  to  validate  the  struc¬ 
ture  of  the  filter’s  error  model  by  tuning  certain  filter  parameters  (noise  strengths  and 
correlation  times)  so  that  the  measurement  residuals  appeared  “correct.”  That  is,  the 
residuals  should  appear  white,  zero-mean,  with  variance  as  predicted  by  the  filter.  If  this 
is  achieved,  and  the  filter’s  error  model  is  sufficiently  valid,  then  the  differences  between 
the  position  and  velocity  estimates  of  the  filter  and  TDAS  should  be  minimized.  The 
justification  for  this  statement  is  that  TDAS-indicated  position  and  velocity  are  derived 
from  precisely  surveyed  points  along  the  length  of  the  track;  therefore,  they  are  best  avail¬ 
able  measurements  of  the  true  quantities.  However  the  poor  geometry  of  the  transponders 
relative  to  the  RSV  trajectory  partially  frustrated  this  goal  by  limiting  the  observability 
of  important  error  states.  The  lack  of  observability  made  the  filter  error  state  estimates 
very  sensitive  to  the  initial  error  state  covariance  estimates.  The  parameters  considered 
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Figure  4.3.  Anomalous  Range  Measurement  Oscillation,  Transponder  211 
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RANGE  RESIDUAL  (feet) 


least  well  known  are:  (1)  the  driving  noise  strengths  the  atmospheric  propagation  delay 
scale  factors,  (2)  the  correlation  coefficients  for  these  driving  noises  for  each  transponder 
pair,  (3)  the  correlation  times  for  these  driving  noises,  (4)  the  driving  noise  strength  for 
the  baro-altimeter  correlated  noise,  and  (5)  the  correlation  time  for  the  baro-altimeter 
correlated  noise. 

Starting  with  the  intial  parameter  values  described  in  Appendix  A,  the  70-state 
filter  was  used  to  process  data  from  three  different  track  runs,  two  south-to-north  and 
one  north-to-south.  Different  combinations  of  tuning  parameters  and  initial  covariance 
estimates  were  tried  for  the  five  items  mentioned  above. 

The  two  sets  of  estimates  sue  compared  at  time  points  corresponding  to  filter  mea¬ 
surement  update  times.  Because  the  TDAS  sampling  times  did  not  coincide  with  the  filter 
measurement  update  times,  the  TDAS  positions  and  velocities  were  interpolated  at  the 
filter  update  times.  This  interpolation  is  accomplished  using  the  MATRIXx  cubic  spline 
function.  After  the  interpolation  is  completed,  the  filter  position  estimates  are  corrected 
for  the  offset  in  the  INS  position  relative  to  the  TDAS  sensor. 

For  the  northward  run,  the  position  correction  was  made  as  follows.  A  positive 
latitude  error  means  that  the  filter  latitude  estimate  is  too  far  north.  A  positive  longitude 
error  means  that  the  filter  longitude  estimate  is  too  far  east.  A  positive  altitude  error 
means  that  the  filter  altitude  estimate  is  too  far  up.  The  position  of  the  INS  relative  to 
the  TDAS  sensor  is  shown  in  Table  4.1.  Since  the  track  heading  relative  to  true  north  is 
constant  throughout  the  length  of  the  track,  these  position  offsets  are  also  constant.  For 
the  northward  and  southward  runs,  the  appropriate  offset  quantities  were  subtracted  from 
the  filter  estimates  in  order  to  translate  the  filter  estimates  to  the  TDAS  sensor  location. 

The  area  where  the  filter  is  most  sensitive  is  the  baro-altimeter  error  model.  The 
separation  of  the  baro-altimeter  error  into  three  separate  error  states  results  in  observ¬ 
ability  problems  when,  as  in  this  case,  the  test  trajectory  is  of  short  duration  and  involves 
only  negligible  altitude  changes.  When  all  three  states  were  included  in  the  filter,  the  sum 
of  the  constant  and  time-correlated  biases  tended  to  grow  increasingly  positive  while  the 
scale  factor  error  grew  increasingly  negative.  The  result  was  the  sum  of  the  three  errors 
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stayed  in  the  range  of  100  to  -100  feet  while  the  individual  error  states  were  divergent. 
This  situation  was  corrected  by  using  only  one  state  to  model  the  baro-altimeter  error 
for  the  ground  test  data.  The  scale  factor  error  state  was  eliminated  and  the  other  two 
error  states  combined  into  a  first-order  Markov  error  state  with  a  correlation  time  of  2000 
seconds  and  a  driving  noise  variance  of  4  ft2.  The  initial  variance  for  this  single  baro- 
altimeter  error  state  and  for  the  INS  altitude  error  state  were  each  reduced  to  2500  ft2. 
Even  with  the  simplified  baro-altimeter  error  model,  the  filter  altitude  estimates  remain 
biased  relative  to  the  TDAS  indicated  altitude,  as  is  shown  in  Figure  4.7. 

The  filter  estimates  are  also  sensitive  to  the  atmospheric  propagation  delay  error 
model  parameters.  The  filter  is  most  sensitive  to  the  assumed  driving  noise  strength  and 
less  sensitive  to  the  assumed  correlation  time  and  driving  noise  correlation  coefficient.  The 
assumed  driving  noise  variance  is  reduced  to  4  ft2  from  the  value  shown  in  Appendix  A. 
This  reduction  had  a  smoothing  effect  on  the  INS  and  transponder  position  error  estimates 
as  well  as  the  atmospheric  delay  estimates. 

The  data  set  discussed  here  is  designated  Track  Rim  E.  Track  Run  E  started  at 
IRIG  time  54826  and  ended  at  IRIG  time  55722.  This  starting  time  is  approximately 
50  minutes  after  the  INS  was  aligned.  Appendix  B  contains  the  plots  of  the  recorded 
data.  The  TDAS  position  and  velocity  data  for  Track  Run  E  is  plotted  in  Figures  B.l 
and  B.2.  The  LN-39  INS  indicated  position,  velocity,  and  acceleration  for  Track  Rim  E 
are  plotted  in  Figures  B.3,  B.4,  and  B.5.  The  baro-altimeter,  transponder  range,  and 
transponder  range-rate  measurements  for  Track  Run  E  are  plotted  in  Figures  B.6  through 
B.12.  The  baro-altimeter  discretization  increment  is  2.5  feet.  The  barometric  altitude  rate 
information  is  synthesized  from  the  barometric  altitude  data. 

The  filter’s  estimates  of  position  and  velocity  are  compared  to  the  TDAS  indicated 
position  and  velocity  in  Figures  4.5  through  4.10.  All  of  these  plots  show  the  filter  estimate 
minus  the  corresponding  TDAS  measurements  The  latitude  difference  varies  between  0.5 
and  -11.5  feet  with  the  larger  magnitudes  in  the  last  half  of  the  run.  The  longitude 
difference  varies  between  16.5  and  -7.0  feet  with  the  larger  magnitudes  in  the  first  half  of 
the  rim.  The  altitude  difference  varies  between  7.8  and  26.5  feet  with  the  larger  magnitudes 
in  the  second  half  of  the  run.  The  midpoint  of  the  run  is  significant  because  that  is  when 
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the  RSV  passed  by  transponder  102,  at  time  55275.  This  is  the  only  transponder  for 
which  the  line-of-sight  changes  significantly  during  the  run.  If  one  assumes  the  TDAS 
measurement  to  be  the  “truth,”  then  the  growth  of  the  latitude  and  altitude  differences 
after  passing  the  midpoint  is  the  opposite  of  what  is  expected  while  the  decrease  in  the 
longitude  difference  is  as  expected.  The  north,  west,  and  up  velocity  differences  are  all 
essentially  zero-mean.  When  the  RSV  is  travelling  at  constant  speed,  the  peak  magnitudes 
of  the  velocity  differences  are  0.6  fps  for  the  north  velocity,  0.3  fps  for  the  west  velocity, 
and  1.7  fps  for  the  up  velocity. 

The  measurement  residuals  are  shown  in  Figures  4.11  through  4.21.  All  residuals 
are  approximately  zero-mean.  The  dashed  lines  on  the  residual  plots  are  the  expected 
residual  ±lcr  bounds  as  computed  by  the  filter.  Only  measurements  within  a  5<r  threshold 
are  accepted  by  the  filter.  For  the  barometric  altitude  and  range-rate  measurements,  the 
residuals  have  a  generally  “white”  appearance  and  most  of  the  residuals  are  within  the 
ler  bounds,  indicating  slightly  over-conservative  filter  tuning  for  these  two  measurement 
types.  For  the  range  measurements,  the  residuals  appear  white  for  some  transponders 
and  less  so  for  others.  The  ratio  of  the  actual  residual  variance  to  the  expected  residual 
variance  also  appears  to  vary  significantly  among  transponders. 

The  filter’s  estimates  of  the  transponder  position  errors  and  the  atmospheric  propa¬ 
gation  delay,  along  with  the  filter  estimated  standard  deviations,  are  shown  in  Figures  4.22 
through  4.31.  The  poor  geometry  of  the  transponders  relative  to  the  track  is  apparent  in 
the  lcr  plots.  There  is  essentially  no  observability  of  the  transponder  vertical  position 
errors.  This  is  because  the  altitude  of  the  RSV  increases  only  30  feet  during  the  run, 
thus  keeping  the  line-of-sight  elevation  angles  for  all  transponders  under  5°.  The  effect 
of  passing  by  transponder  102  is  apparent  in  most  of  the  standard  deviation  plots;  this  is 
where  the  l<r  values  decrease  most  rapidly. 

The  magnitudes  of  the  estimated  position  errors  remain  resonable  values  throughout 
the  run,  remaining  under  7.5  feet.  Some  of  the  position  error  estimates  appear  to  converge 
on  a  relatively  steady  value  in  the  second  half  of  the  run  while  others  continue  to  vary. 
Of  particular  interest  is  the  spike  in  the  estimated  west  position  error  for  transponder 
102.  This  spike  seems  to  be  correlated  with  spikes  in  the  estimated  position  errors  of  the 
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other  transponders,  particularly  transponders  005  and  212.  These  rapid  variations  indicate 
the  filter’s  redistribution  of  the  measurement  residual  error  during  the  period  where  it  is 
gaining  the  most  information  about  the  actual  position  errors  for  transponder  102. 

The  magnitudes  of  the  estimated  atmospheric  propagation  delays,  after  seeding  for 
range,  correspond  to  linear  range  errors  between  2.5  feet  and  -2.5  feet.  Although  the 
true  value  of  the  propagation  delay  must  always  be  positive,  the  negative  estimates  are 
possible  for  two  reasons.  The  first  is  the  inability  of  the  filter  to  distinguish  between  the 
atmospheric  delay  and  the  interrogator  range  calibration  bias  within  the  short  duration  of 
the  run.  The  second  is  the  possibility  of  error  in  the  open  loop  correction  for  propagation 
delay  made  by  the  CIRIS  interrogator.  The  1  <t  plots  for  the  atmospheric  delay  error  states 
appear  to  be  converging  to  steady  values. 
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Figure  4.6.  Longitude  Difference,  Track  Run  E 
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Figure  4.7.  Altitude  Difference,  Track  Run  E 
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Figure  4.8.  North  Velocity  Difference,  Track  Run  E 
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Figure  4.9.  West  Velocity  Difference,  Track  Run  E 
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Figure  4.10.  Vertical  Velocity  Difference,  Track  Run  E 
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Figure  4.11.  Baro-Altitude  Residuals,  Track  Run  E 
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Figure  4.12.  Transponder  005  Range  Residuals,  Track  Run  E 
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Figure  4.13.  Transponder  005  Range-Rate  Residuals,  Track  Run  E 
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Figure  4.14.  Transponder  102  Range  Residuals,  Track  Run  E 
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Figure  4.15.  Transponder  102  Range-Rate  Residuals,  Track  Run  E 
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Figure  4.16.  Transponder  181  Range  Residuals,  Track  Run  E 
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Figure  4.17.  Transponder  181  Range-Rate  Residuals,  Track  Run  E 
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Figure  4.19.  Transponder  212  Range-Rate  Residuals,  Track  Run  E 
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Figure  4.22.  Estimated  Position  Error,  Trans  005,  Track  Run  E 
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Figure  4.23.  Estimated  Position  Std  Dev,  Trans  005,  Track  Run  E 
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Figure  4.25.  Estimated  Position  Std  Dev,  Trans  102,  Track  Run  E 
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Figure  4.26.  Estimated  Position  Error,  Trans  181,  Track  Run  E 
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Figure  4.27.  Estimated  Position  Std  Dev,  Trans  181,  Track  Rim  E 
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Figure  4.29.  Estimated  Position  Std  Dev,  Trans  212,  Track  Run  E 
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4-3  Summary 

This  chapter  presents  the  results  of  using  the  Advanced  CIRIS  filter  to  process  real 
CIRIS  data  (measurements  and  INS  output)  from  a  slow-speed  ground  test.  The  objective 
is  to  validate  the  filter’s  error  model  by  analysis  of  the  measurement  residual  characteris¬ 
tics  and  by  comparison  of  the  filter’s  position/velocity  estimates  with  those  of  the  track 
data  acquisition  system.  Range  measurements  from  one  of  the  six  transponders  used  in 
the  test  contain  an  anomalous  oscillation  (90  second  period),  so  measurements  from  this 
transponder  are  not  used.  In  general,  the  baro-altimeter  and  range-rate  residuals  ap¬ 
pear  zero-mean  and  white,  with  the  expected  residual  variance  achievable  by  tuning  the 
measurement  noise  parameters.  The  range-residuals,  although  approximately  zero-mean, 
exhibit  short  duration  biases,  and  the  actual  residual  variance  appears  to  vary  significantly 
among  transponders.  The  errors  in  the  filter’s  position  estimates,  relative  to  the  TDAS 
measurements,  fall  mostly  in  the  ±20  ft  range  and  exhibit  biases.  The  corresponding  veloc¬ 
ity  errors  are  mostly  zero-mean  with  peak  magnitudes  in  the  ±0.3  fps  range,  although  the 
up-velocity  error  magnitudes  are  significantly  larger.  The  problem  with  the  up-velocity 
magnitudes  is  most  likely  due  to  the  vibration  environment.  The  filter’s  estimates  of  the 
transponder  survey  position  errors  and  atmospheric  propagation  scale  factor  appear  rea¬ 
sonable,  although  the  associated  1<t  plots  provide  definite  indications  of  near-zero  vertical 
position  error  observability.  The  overall  results  are  interpreted  as  indicating  the  current 
filter  structure  is  reasonable,  but  the  adequacy  of  the  measurement  models  requires  further 
investigation. 
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V.  70-State  Filter  Evaluation:  Flight  Test  Data 

In  this  chapter,  the  performance  of  the  70-state  filter  is  evaluated  using  data  collected 
during  a  C-141  aircraft  flight  over  the  Yuma  Test  Range. 

5.1  Data  Collection  Procedure 

The  data  used  in  this  analysis  was  collected  on  12  September  1989.  The  CIRIS 
LN-39  INS,  transponder  interrogator,  baro-altimeter,  and  HP-1000  computer  were  placed 
on  a  C-141  aircraft  that  flew  from  Holloman  AFB  to  Yuma,  Arizona  and  back.  During 
the  portion  of  the  flight  over  the  Yuma  Test  Range,  the  aircraft  was  tracked  by  a  laser 
ranging  system  (LRS),  from  which  aircraft  position  and  velocity  were  determined.  The 
LRS  measurements  are  known  to  be  about  twice  as  accurate,  on  average,  than  CIRIS 
I  filter  estimates  [7].  The  aircraft  trajectory,  as  determined  by  the  LRS,  is  plotted  in 
Figure  5.1.  Appendix  C  contains  more  detailed  plots  of  both  the  LRS  and  LN-39  INS 
indicated  position  and  velocity  for  the  Yuma  Test  Range  portion  of  the  flight. 

During  the  entire  flight,  the  CIRIS  equipment  recorded  measurements  from  twenty- 
two  different  transponders.  Since  the  70-state  filter  software  can  accommodate  only  20 
transponders,  measurements  from  the  two  least  used  transponders  are  not  processed  by 
this  filter.  The  remaining  twenty  transponders  are  paired  up  such  that  the  state  switching 
procedure  described  in  Chapter  III  is  needed  only  once  during  the  Yuma  Test  Range 
portion  of  the  flight. 

5.2  70-State  Filter  Performance 

For  this  flight  test,  the  LRS  indicated  position  and  velocity  is  used  as  the  “truth” 
by  which  the  the  accuracy  of  the  70-state  filter  is  evaluated.  The  initial  state  covariance 
estimates  and  other  tuning  parameters  used  in  this  test  are  exactly  those  described  in 
Appendix  A.  The  recording  of  LN-39  INS  data  and  transponder  and  baro-altimeter  mea¬ 
surements  began  prior  to  take-off  and  continued  until  the  aircraft  landed.  The  70-state 
filter  is  used  to  process  this  data  from  just  prior  to  take-off  until  the  aircraft  leaves  the 
Yuma  area.  This  requires  almost  11  hours  of  CPU  time  on  a  MicroVAX  III  computer  to 
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process  3.6  hours  of  real-time  data.  This  figure  includes  the  large  amount  of  I/O  required 
to  read  the  input  data  and  save  calculated  results  for  plotting. 

The  differences  between  the  filter  estimated  position  and  velocity  and  the  LRS  indi¬ 
cated  position  and  velocity  are  plotted  in  Figures  5.2  through  5.7.  Compensation  is  made 
for  the  position  offset  between  the  LN-39  INS  and  the  laser  retro-reflector. 

The  baro-altimeter  measurement  residuals  are  plotted  in  Figure  5.8.  The  transpon¬ 
der  range  and  range-rate  residuals  for  the  three  most  frequently  interrogated  transponders, 
160,  163,  and  226,  are  shown  in  Figures  5.9  through  5.14.  The  bounding  lines  on  these 
plots  are  the  residual  ±l<r  values  predicted  by  the  filter.  The  baro-altimeter  residuals  are 
near  zero-mean  and  have  a  true  variance  much  smaller  than  the  expected  variance  except 
when  the  aircraft  is  climbing  or  diving.  During  altitude  changes,  the  residual  variance 
increases  dramatically.  These  plots  indicate  the  filter  baro-altimeter  measurement  noise 
parameter  is  set  too  high.  It  should  probably  be  adjusted  as  a  function  of  the  vertical 
velocity.  The  range  measurement  for  all  three  transponders  look  nearly  ideal.  They  are 
near  zero-mean,  appear  white  most  of  the  time,  and  have  actual  variances  close  to  the 
expected  variances.  The  range-rate  residuals  from  all  three  transponders  do  not  possess 
these  ideal  characteristics.  They  appear  near  zero-mean  overall,  but  they  exhibit  time 
correlated  variations  and  have  actual  variances  much  larger  than  the  expected  variances 
(some  large  magnitude  residuals  are  not  shown  on  the  plot). 

The  filter  residual  tolerance  was  set  to  5 <j.  During  the  portion  of  the  flight  from  Hol¬ 
loman  to  Yuma  (IRIG  times  48660  to  53000)  the  filter  received  2,946  barometric  altitude 
measurements,  2,946  transponder  range  measurements,  and  2,946  transponder  range-rate 
measurements.  Of  these,  none  of  the  barometric-altitude  measurements  were  rejected,  97 
(3.293%)  of  the  range  measurements  were  rejected,  and  135  (4.582%)  of  the  range-rate 
measurements  were  rejected.  Fifty-two  of  the  rejected  range  measurements  were  from 
transponder  218.  During  the  portion  of  the  flight  over  the  Yuma  Test  Range  (IRIG  times 
53000  to  63000)  the  filter  received  5,663  barometric  altitude  measurements,  5,663  transpon¬ 
der  range  measurements,  and  5,663  transponder  range-rate  measurements.  Of  these,  one 
(0.018%)  of  the  barometric-altitude  measurements  was  rejected,  83  (1.466%)  of  the  range 
measurements  were  rejected,  and  1,656  (29.242%)  of  the  range-rate  measurements  were  re- 


5-2 


jected.  Forty-two  of  the  rejected  range  measurements  were  from  transponder  218.  However 
the  majority  of  the  range  measurements  from  this  particular  transponder  were  accepted 
by  the  filter  in  both  portions  of  the  flight. 

The  rejection  of  the  range-rate  measurements  is  correlated  with  the  aircraft  turn¬ 
ing  maneuvers.  This  may  indicate  that  modeling  the  range-rate  measurement  process  as 
occurring  instantaneously  is  an  over-optimistic  assumption.  The  way  the  range-rate  mea¬ 
surements  are  made  by  the  Cubic  RRS  actually  involves  calculating  a  delta-range  (distance 
units)  over  a  short  (approximately  0.5  seconds)  measurement  time  interval.  The  delta- 
range,  after  dividing  by  the  measurement  time  interval,  represents  the  average  range-rate 
during  the  interval.  During  aircraft  turning  maneuvers,  because  the  INS  is  not  at  the 
aircraft’s  center  of  rotation,  this  average  range-rate  may  not  accurately  reflect  the  true 
instantaneous  range-rate  at  the  end  of  the  measurement  interval  when  the  measurement 
is  time-tagged.  This  possibility  requires  further  research;  it  may  necessitate  a  redesign  of 
the  filter’s  range-rate  measurement  model. 

5.3  CIRIS  II  Filter  Performance 

The  CIRIS  II  filter  position  and  velocity  estimates  are  available  for  this  same  flight, 
so  it  is  useful  to  compare  these  estimates  to  the  LRS  indicated  position  and  velocity.  This 
permits  an  evaluation  of  the  performance  of  Advanced  CIRIS  relative  to  CIRIS  II.  The 
CIRIS  II  estimates  used  to  make  the  plots  in  this  section  are  not  corrected  for  the  INS 
to  laser-retroreflector  offset,  however  this  type  of  correction  made  minimal  difference  in 
the  corresponding  Advanced  CIRIS  filter  plots.  The  CIRIS  II  transponder  measurement 
residuals  are  also  available  so  they  may  be  compared  to  the  corresponding  Advanced  CIRIS 
residuals.  The  CIRIS  II  filter  does  not  use  separate  barometric  altitude  measurements. 

The  CIRIS  II  latitude,  longitude,  and  altitude  errors,  relative  to  the  LRS  measure¬ 
ments,  are  shown  in  Figures  5.15,  5.16,  and  5.17.  These  plots  exhibit  spikes  similar  to  the 
ones  in  the  Advanced  CIRIS  plots,  but  the  average  error  magnitudes  are  larger,  especially 
for  the  altitude  error.  The  CIRIS  II  north,  west,  and  up  velocity  errors,  relative  to  the 
LRS  measurements,  are  shown  in  Figures  5.18,  5.19,  and  5.20.  These  plots  also  contain 
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spikes  similar  to  the  ones  in  the  Advanced  CIRIS  plots  and  the  average  error  magnitudes 
are  roughly  the  same  as  those  of  Advanced  CIRIS. 

The  CIRIS  II  range  and  range-rate  measurement  residuals  for  transponders  160, 163, 
and  226  are  shown  in  Figures  5.21  through  Figures  5.26.  Compared  to  the  the  correspond¬ 
ing  Advanced  CIRIS  range  residuals,  the  CIRIS  II  range  residuals  are  obviously  biased  and 
have  a  larger  variance.  The  CIRIS  II  range-rate  residuals  are  also  biased,  but  they  have 
an  obviously  smaller  variance  than  the  Advanced  CIRIS  range-rate  residuals. 
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Figure  5.2.  Latitude  Difference,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.5.  North  Velocity  Difference,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.7.  Up  Velocity  Difference,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.8.  B mo- Altitude  Meas.  Residuals,  Yuma  Flight,  Adv.  CIRJS 
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Figure  5.9.  Trans.  160  Range  Meas.  Residuals,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.10.  Trans.  163  Range  Meas.  Residuals,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.11.  Trans.  226  Range  Meas.  Residuals,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.12.  Trans.  160  Range-Rate  Meas.  Residuals,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.13.  Trans.  163  Range-Rate  Meas.  Residuals,  Yuma  Flight,  Adv.  CIRIS 
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Figure  5.14.  Trans.  226  Range-Rate  Meas.  Residuals,  Yuma  Flight,  Adv.  CIRJS 
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Figure  5.15.  Latitude  Difference,  Yuma  Flight,  CIRIS  II 
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Figure  5.16.  Longitude  Difference,  Yuma  Flight,  CIRIS  II 
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Figure  5.17.  Altitude  Difference,  Yuma  Flight,  CIRIS  II 
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Figure  5.18.  North  Velocity  Difference,  Yuma  Flight,  CIRIS  II 
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Figure  5.19.  West  Velocity  Difference,  Yuma  Flight,  CIRIS  II 
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Figure  5.20.  Up  Velocity  Difference,  Yuma  Flight,  CIRIS  II 
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Figure  5.21.  Transponder  160  Range  Measurement  Residuals,  CIRIS  II 
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Figure  5.22.  Transponder  163  Range  Measurement  Residuals,  CIRIS  II 
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Figure  5.23.  Transponder  226  Range  Meas.  Residuals,  CIRIS  II 
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Figure  5.24.  Transponder  160  Range-Rate  Meas.  Residuals,  CIRIS  II 
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Figure  5.25.  Transponder  163  Range-Rate  Meas.  Residuals,  CIRIS  II 
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Figure  5.26.  Transponder  226  Range-Rate  Meas.  Residuals,  CIRIS  II 
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5-4  Summary 

This  chapter  presents  the  results  of  using  the  Advanced  CIRIS  filter  to  process  real 
CIRIS  data  recorded  during  a  C-141  aircraft  flight  over  the  Yuma  Test  Range.  Dur¬ 
ing  flight,  the  aircraft  was  tracked  by  an  independent  laser  ranging  system  (LRS)  whose 
measurements  are  known  to  be  roughly  twice  as  accurate,  on  average,  as  CIRIS  I  filter  esti¬ 
mates.  Use  of  this  data  provides  additional  filter  model  validity  information  through  anal¬ 
ysis  of  residual  characteristics.  In  general,  the  Advanced  CIRIS  filter’s  measurement  resid¬ 
uals  appear  reasonably  white  and  zero-mean.  The  actual  variance  of  the  baro-altimeter 
measurement  residuals  increases  significantly  during  aircraft  altitude  changes.  The  actual 
variance  of  the  range-rate  measurement  residuals  appears  correlated  with  aircraft  ma¬ 
neuvers;  the  large  residual  magnitudes  during  turning  maneuvers  result  in  an  excessive 
proportion  of  the  range-rate  measurements  being  rejected  by  the  filter.  By  comparison, 
the  CIRIS  II  range  measurement  residuals  exhibit  a  significantly  larger  true  variance  than 
those  of  Advanced  CIRIS,  and  they  also  exhibit  long  duration  biases.  However  the  CIRIS 
II  range-rate  measurement  residuals,  while  also  biased,  have  an  obviously  smaller  actual 
variance  than  those  of  Advanced  CIRIS. 

The  availability  of  CIRIS  II  filter  estimates  for  this  same  flight  permits  am  estimation 
accuracy  comparison  between  Advanced  CIRIS  and  CIRIS  II,  both  relative  to  the  LRS 
measurements.  On  average,  the  relative  errors  in  the  Advanced  CIRIS  position  and  velocity 
estimates  are  smaller  than  those  of  CIRIS  II.  The  advantage  of  Advanced  CIRIS  is  most 
noticeable  in  the  altitude  error  plots. 
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VI.  Proposed  Fixed-Interval  Smoother 


This  chapter  describes  a  fixed-interval  smoother  algorithm  whose  simple  structure 
makes  it  a  good  companion  program  for  the  current  MSOFE-based  Advanced  CIRIS  filter 
software.  The  smoother  software  was  not  implemented  in  this  thesis  effort  due  to  lack  of 
•  time  and  concern  over  data  storage  requirements. 

6.1  Fixed- Interval  Smoothing 

The  term  smoother  is  used  to  describe  a  type  of  estimator  that  makes  use  of  some 
or  all  measurements  taken  after  the  time  of  interest  as  well  as  all  of  the  measurements 
preceeding  that  time  [9].  This  characteristic  prevents  the  real-time  computation  of  the 
smoothed  estimates,  but  it  provides  the  potential  for  more  accurate  estimates.  Given  a 
time  interval  of  interest,  [ta,tf\  and  a  time  ti  in  this  interval,  a  Kalman  filter  propagates  the 
conditional  probability  density  function  shown  in  Equation  (2.24)  to  yield  the  expectation 
of  x(fj)  conditioned  on  Z;,  the  measurement  history  through  time  U.  This  is  expressed  as 


x(tf)  = 

£{x(ti)  |  Z(tj,u>fe)  =  Zi} 

(6.1) 

x(t“)  = 

E{x(ti)  I  Z(t;_i  ,U>k)  =  Z;_i} 

(6.2) 

In  contrast,  the  smoothed  estimate  at  time  ti  is  conditioned  on  the  measurement  history 
through  time  tj,  where  >  f  i.  This  is  written  as 

xiU/tj)  =  E{x{ti)  |  Z(tj,u>fc)  =  Z;},  j  >  i  (6.3) 

A  fixed-interval  smoother  uses  all  measurements  through  the  final  time.  In  other  words, 
tj  =  tf.  The  other  types  of  smoothers  are  fixed-point  and  fixed-lag.  This  chapter  is  not 
intended  to  present  the  analytic  development  of  smoother  equations.  For  those  details,  the 
reader  should  consult  references  [9,10]. 

A  fixed-interval  smoothing  algorithm  may  be  implemented  as  a  two-pass  Kalman 
filtering  procedure:  a  standard  (or  extended)  Kalman  filter  for  the  forward  pass  and  an 
inverse-covariance  (or  extended  inverse-covariance)  Kalman  filter  for  the  backward  pass. 
A  smoothing  algorithm  may  be  further  characterized  as  discrete-time  or  continuous-time, 
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depending  on  the  underlying  state  dynamics  propagation  model  (state  transition  matrix  or 
matrix  differential  equation,  respectively).  The  MSOFE-based  Advanced  CIRIS  Kalman 
filter  uses  a  continuous  time  model  so,  ideally,  a  companion  smoothing  program  should  use  a 
continuous-time  smoothing  algorithm.  However,  the  continuous-time  smoother  equations 
are  quite  complicated  and  difficult  to  implement  [10].  For  this  reason,  the  simple  formula¬ 
tion  for  a  discrete-time  fixed-interval  smoother  developed  by  Meditch  [10]  is  adapted  for 
use  with  the  Advanced  CIRIS  Kalman  filter.  This  algorithm  is  also  described  in  reference 
[9:14-15]. 

The  smoothing  procedure  is  as  follows.  During  the  forward  filter  pass,  the  quantities 
x(t“),  P (£“),  x(tf ),  and  P (tf )  are  saved  for  all  measurement  times  in  the  interval  [t0,  t/]. 
The  dynamics  matrix  F[x(t*)]  must  also  be  saved  at  these  times  for  use  in  computing  the 
state  transition  matrix.  Starting  from  the  boundary  condition 

=  x(t£)  (6.4) 

the  smoothed  estimate  is  generated  backward  in  time  using 

x(ti/tf)  =  x{tf)  +  A  (ti)  [x(ti+1/t/)  -  x(t-+1)]  (6.5) 

where  the  “smoothing  estimator  gain  matrix”  A  (ti)  is  given  by 

A  (ti)  =  P(tl+)#T(ti+i,t,)P~1(fr+i)  (6-6) 

The  state  transition  matrix  #(t,+i,t,),  whose  transpose  is  required  in  Equation  (6.6),  is 
the  solution  to  the  differential  equation 

#(*,*<_!)  =  F[x(t/<i_i)]#(Mi_i),  t€[ti-i,ii)  (6.7) 

If  the  time  interval  [t;_i,t,]  is  sufficiently  small,  so  that  F[x(t/fi_i)]  is  essentially  constant 
over  the  interval,  then  #(t,+i,f,)  may  be  calculated  using 

=  exp  [F[x(t,++1)](ti+i  -  ti)}  (6.8) 

Analogous  to  the  forward-pass  filter  covariance  matrix,  the  covariance  matrix  for  the 
smoothed  estimates  is  P(ti/£/),  This  smoother  covariance  matrix  retains  the  interpreta¬ 
tion  as  the  covariance  of  the  zero-mean  Gaussian  estimation  error  [x(t,/t/)  -  x(tj)].  The 
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smoother  covariance  is  generated  backwards  starting  from  the  boundary  condition 

P  (U/tf)  =  P  (t+)  (6.9) 

using  the  recursion 

P (ti/tf)  =  P (tf)  +  A(tt)  [P(ti+1/t/)  -  P(t-+1)]  At(U)  (6.10) 

6.2  Application  to  the  Advanced  CIRIS  Filter 

The  smoothing  procedure  described  in  the  previous  section  may  require  large  amounts 
of  disk  space,  depending  on  the  state  vector  length  and  the  number  of  measurement  times 
in  the  fixed-interval.  Consider  a  one  hour  interval  with  a  measurement  vector  processed 
every  second;  the  interval  contains  3600  measurement  times.  The  quantities  x(t~ ),  P (t~), 
x(tf),  P (tf)  and  F[x(t+)]  must  be  saved  at  each  measurement  time.  Because  the  P 
matrix  is  symmetric,  only  the  upper  triangular  part  need  be  stored.  For  a  state  vector 
length  of  70,  x(f,~)  and  x(tf)  each  require  70  64-bit  data  words,  the  upper  triangles  of 
P(tf)  and  P(#t+ )  each  require  2,485  64-bit  data  words,  and  F[x(t,+ )]  requires  4,900  64- 
bit  data  words.  The  result  is  10,010  data  words  (80,080  bytes)  requiring  storage  at  each 
measurement  time.  The  total  (unformatted)  disk  storage  required  for  the  one  hour  interval 
is  over  288  megabytes.  Storage  of  time-tag  information  requires  another  28,800  bytes. 
While  this  amount  of  disk  space  may  be  available  on  a  dedicated  computer,  such  resources 
are  not  yet  generally  available. 

Storage  requirements  can  be  significantly  reduced  by  taking  advantage  of  the  special 
structure  of  the  F[x(t+ )]  matrix  in  the  Advanced  CIRIS  filter.  As  indicated  in  Appendix  A, 
Tables  A. 7  through  A. 9,  the  the  majority  of  the  F  matrix  elements  are  constant,  most 
of  them  zero.  Since  all  of  the  time- varying  elements  are  contained  in  the  upper  left 
28-by-28  submatrix,  only  this  submatrix  should  be  stored  at  each  measurement  time.1 
The  remaining  42,  constant,  non-zero,  diagonal  elements  need  be  saved  only  once.  This 
alternative  procedure  reduces  the  required  F  matrix  storage  at  each  measurement  time  to 
784  64-bit  data  words  and  reduces  the  total  (unformatted)  storage  requirement  for  the  one 
hour  interval  to  170  megabytes. 

‘For  even  greater  reduction*,  only  the  time  varying  element*  should  be  stored.  This  would  require  a 
slight  increase  in  the  complexity  of  the  storage  procedure. 
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The  other  costly  aspect  of  the  proposed  smoother  algorithm  is  computation  time.  In 
addition  to  the  normal  matrix  arithmetic,  the  procedure  requires  a  matrix  inversion  and 
a  matrix  exponential  for  each  measurement  time.  A  search  of  available  references  did  not 
reveal  any  special  procedures  for  finding  the  symmetric  inverse  of  a  symmetric  matrix  using 
only  the  upper  triangle  of  the  matrix.  Numerical  procedures  exist  for  computing  matrix 
exponentials  [5],  but  they  must  be  carefully  implemented  to  ensure  accurate  results.  In  the 
case  of  Advanced  CIRIS,  the  block  diagonal  structure  of  the  F  matrix  simplifies  calculation 
of  the  state  transition  matrix.  Because  the  lower  right  42-by-42  submatrix  is  constant  and 
diagonal,  the  exponential  of  this  submatrix  is  simply  another  constant,  diagonal,  submatrix 
formed  by  exponentiating  each  of  the  original  diagonal  elements.  The  exponential  of  the 
upper  left  28-by-28  submatrix  may  be  calculated  using  the  algorithm  described  in  reference 
[5:396-400].  The  remaining  (upper  right  and  lower  left)  submatrices  are  identically  zero 
so  their  exponentials  are  also  zero.  Care  must  be  taken  to  ensure  numerically  accurate 
results  at  each  step  in  the  backward  recursion.  Since  the  non-recursive  calculation  of  the 
smoothing  estimator  gain  matrix  depends  only  on  the  previously  recorded  data  for  each 
measurement  time,  this  calculation  is  not  subject  to  the  cumulative  effects  of  round-off 
error.  However,  the  recursive  calculations  of  the  smoothed  state  estimates  and  smoothed 
covariance  matrices  are  subject  to  such  cumulative  effects. 

The  adaptation  of  a  discrete-time  smoothing  algorithm  to  a  continuous-time  problem 
contains  the  implicit  assumption  that  the  dynamic  matrix  F  is  relatively  constant  over 
the  time  interval  between  measurements.  If  this  assumption  is  not  true,  then  the  state 
transition  matrix  calculated  at  each  measurement  time  will  fail  to  represent  the  actual 
state  dynamics  between  measurement  times.  With  typical  CIRIS  data,  the  measurements 
occur  about  one  second  apart,  but  there  is  no  absolute  upper  bound  on  this.  Since  the  time 
varying  eigenvalues  of  the  Advanced  CIRIS  dynamics  matrix  all  have  time  constants  greater 
than  about  42  minutes  (half  the  Schuler  period),  and  the  rate  of  change  of  these  eigenvalues 
is  relatively  slow  for  cargo  type  aircraft  [7],  the  assumption  that  F  is  essentially  constant 
seems  resonable.  Another  area  of  concern  is  whether  or  not  the  proposed  fixed-interval 
smoothing  algorithm  is  compatible  with  the  ad  hoc  procedure  described  in  Chapter  III  for 
the  time-sharing  of  transponder  error  states.  It  is  not  readily  apparent  that  the  backward 
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recursion  can  proceed  correctly  when  the  forward  pass  covariance  matrix  is  discontinuous  at 
times  t~ .  Until  this  is  resolved,  the  smoothing  algorithm  should  not  be  used  in  conjunction 
with  time-shared  transponder  error  states. 

6.3  Summary 

This  chapter  introduces  the  concept  of  a  fixed-interval  smoother  as  an  optimal  es¬ 
timator  based  on  a  separate  forward-pass  and  backward-pass  filters.  An  adaptation  of 
simple  discrete-time  dynamics  formulation  for  the  backward-pass  filter  is  proposed  for 
use  with  the  continuous-time  dynamics  model  in  the  forward-pass  Advanced  CIRIS  filter. 
The  reason  for  this  is  to  avoid  the  computational  complexity  of  a  true  continuous-time 
dynamics  backward-pass  filter.2  Because  the  backward  pass  requires  use  of  large  amounts 
of  data  stored  during  the  forward  pass,  use  of  the  smoother  for  very  long  data  sets  may 
not  be  practical  due  to  disk  storage  limitations.  The  proposed  smoother  was  not  imple¬ 
mented  in  this  thesis  effort  due  to  lack  of  time,  but  it  holds  promise  as  a  means  of  achieving 
additional  estimation  accuracy  beyond  that  of  the  Advanced  CIRIS  filter  alone. 


JMSOFE  itself  does  not  provide  the  subroutines  necessary  to  implement  backward  pass  filters. 
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VII.  Conclusions  and  Recommendations 

This  chapter  restates  the  major  goals  of  this  thesis  effort,  evaluates  progress  made 
towards  these  goals,  and  provides  recommendations  for  further  research  regarding  the 
Advanced  CIRIS  filter. 

7. 1  Conclusions 

The  first  goal  is  to  improve  the  structure  and  utility  of  the  MSOFE-based  software 
which  implements  the  filter  described  in  this  thesis.  Although  the  preceding  chapters  do 
not  specifically  address  the  filter  software  in  detail,  it  is  the  software  which  consumed  the 
majority  of  the  effort  required  for  this  research.  The  software  developed  by  Solomon  [13] 
was  procedurally  correct  but  contained  numerous  embedded  constants  unique  to  the  CIRIS 
data  set  used  in  his  research.  This  software  was  completely  rewritten  to  eliminate  most 
embedded  constants,  allow  tuning  parameters  to  be  changed  without  recompiling,  increase 
modularity,  and  greatly  simplify  control  structures. 

The  second  and  primary  goal  is  to  evaluate,  using  read  CIRIS  data,  the  performance 
of  a  refined  version  of  the  70-state  filter  developed  by  Solomon.  Refinements  to  the  filter 
include  modeling  the  range  measurement  errors  due  to  atmospheric  propagation  delays 
as  being  range  dependent,  modeling  the  cross-correlation  of  the  atmospheric  error  driv¬ 
ing  noises  for  different  transponders,  and  using  barometric-altitude  measurements.  The 
development  of  the  range  and  range-rate  measurement  error  models  is  reaccomplished  to 
illustrate  the  procedure  for  compensating  for  the  offset  between  the  interrogator  antenna 
and  the  INS.  The  filter  performance  was  evaluated  using  two  separate  data  sets. 

The  first  data  set  came  from  a  series  of  slow  speed  ground  tests.  The  reason  for  using 
this  data  is  the  availability  of  independent  position  and  velocity  measurements  provided 
by  the  track  data  acquisition  system  (TDAS).  In  addition  to  evaluating  model  validity 
by  analysis  of  residual  characteristics,  these  independent  measurements  are  considered 
accurate  enough  to  justify  the  adjustment  of  filter  tuning  parameters  in  order  to  minimize 
the  differences  between  the  filter  position/ velocity  estimates  and  the  TDAS  measurements. 
However,  the  poor  geometry  of  the  usable  transponders  relative  to  the  vehicle  trajectory 
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limited  the  observability  of  most  filter  error  states.  This  was  apparent  by  the  filter- 
calculated  error  state  variances  remaining  relatively  large  throughout  the  test.  Use  of 
this  data  did  indicate  that  modeling  the  baro-altitude  error  with  three  separate  states, 
as  is  suggested  in  the  Litton  LN-39  Systems  Engineering  Report,  is  overly  complex  for 
nearly  constant  altitude  trajectories.  For  the  particular  track  run  discussed  in  Chapter  IV, 
the  peak  position  estimation  errors,  relative  to  the  track  measurements,  were  12  feet  in 
latitude,  17  feet  in  longitude,  and  29  feet  in  altitude.  In  retrospect,  the  time  spent  trying 
to  tune  the  filter  using  this  data  set  was  excessive  considering  the  observability  problems. 
However,  use  of  this  data  did  provide  insights  regarding  the  validity  of  the  error  model  on 
which  the  filter  is  based. 

The  second  data  set  came  from  a  C-141  flight.  This  Sight  involved  many  maneu¬ 
vers  at  different  altitudes  so,  except  for  certain  transponder  position  errors,  there  are  no 
observability  problems  inherent  in  the  recorded  measurements.  The  reason  for  using  this 
data  is  the  availability  of  independent  position  and  velocity  measurements  provided  by  a 
laser  ranging  system  (LRS).  The  results  from  this  data  set  more  are  difficult  to  interpret. 
For  the  Advanced  CIRIS  filter,  the  baro-altimeter  and  transponder  range  measurement 
residuals  appear  “white,”  zero-mean,  with  true  variance  near  the  filter-predicted  variance. 
The  transponder  range-rate  measurement  residuals  are  near  zero-mean,  but  they  exhibit 
some  time  correlation  and  their  actual  variance  is  larger  than  the  filter-predicted  variance. 
The  large  residuals  are  likely  correlated  with  aircraft  turns,  but  more  analysis  is  needed  to 
prove  this.  The  most  likely  explanation  is  that  the  range-rate  measurements  are  not  made 
instantaneously;  they  represent  an  average  range-rate  over  a  short  (  1  second  or  less)  time 
interval.  However,  the  “Advanced”  CIRIS  filter  models  them  as  instantaneous. 

The  comparison  of  the  Advanced  CIRIS  filter  position  and  velocity  estimates  to 
the  LRS  measurements  provides  encouraging  results  but  also  raises  questions  about  the 
comparison  process.  The  filter  latitude  estimate,  though  biased,  differs  from  the  laser  data 
less  than  15  feet  during  most  of  the  flight.  The  filter  longitude  estimate,  also  biased  and 
quite  noisy,  differs  from  the  laser  data  less  than  20  feet  during  most  of  the  flight.  The 
filter  altitude  estimate,  again  biased  and  noisy,  differs  from  the  laser  data  less  than  35 
feet  during  most  of  the  flight.  The  filter  north,  west,  and  up  velocity  estimation  errors 
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are  Jill  near  zero-mean  with  most  errors  in  the  range  +3  to  -3  fps.  However  fill  six  plots, 
especially  the  velocity  difference  plots,  exhibit  large  spikes. 

The  laser  measurement  data  file  provides  latitude,  longitude,  altitude,  and  the  three 
velocities  at  each  laser  measurement  time.  The  measurements  at  each  time  are  provided  by 
the  single  laser  instrument  closest  to  the  aircraft  [7].  The  LRS  data  file  is  thus  a  sequence  of 
groups  of  position  and  velocity  measurements  provided  by  different  laser  instruments.  Since 
there  is  no  averaging  or  filtering  of  the  LRS  data,  the  position  and  velocity  measurements 
may  change  by  relatively  large  amounts  when  the  laser  instrument  is  changed.  In  order  to 
make  the  position  and  velocity  comparisons,  the  laser  measurements  must  be  interpolated 
at  the  filter  measurement  times.  The  use  of  cubic  splines  to  interpolate  discontinuous  data 
frequently  results  in  erroneous  cubic  (polynomial)  oscillations  at  the  points  of  discontinuity. 
If  the  interpolated  data  is  resampled  at  an  oscillation  point,  the  result  is  erroneous.  This 
is  a  possible  cause  of  the  spikes  apparent  in  the  difference  plots. 

The  CIRIS  II  filter  position  estimates,  velocity  estimates,  and  transponder  measure¬ 
ment  residuals  were  available  for  this  flight.  This  CIRIS  II  data  was  compared  to  the  LRS 
measurements  to  gain  insight  into  the  relative  performance  of  CIRIS  II  and  Advanced 
CIRIS.  The  CIRIS  II  range  residuals  exhibit  constant  biases  and  significantly  larger  vari¬ 
ance  than  those  of  the  Advanced  CIRIS  filter.  The  CIRIS  II  range-rate  residuals  exhibit 
constant  biases  but  they  have  a  much  smaller  variance  them  those  of  the  Advanced  CIRIS 
filter.  The  presence  of  the  transponder  position  error  states  in  the  Advanced  CIRIS  filter 
is  the  probable  reason  for  the  absence  of  biases  in  its  residuals.  The  problem  with  the 
large  variance  of  the  Advanced  CIRIS  range-rate  measurement  residuals  appears  to  be  a 
modeling  problem  rather  than  a  parameter  tuning  problem. 

The  comparison  of  the  CIRIS  II  position  and  velocity  estimates  to  the  LRS  measure¬ 
ments  shows  the  CIRIS  II  position  errors  to  be  larger  than  those  of  Advanced  CIRIS  while 
the  velocity  errors  appear  very  similar.  The  CIRIS  II/LRS  difference  plots  exhibit  spikes 
similar  to  the  ones  in  the  Advanced  CIRIS/LRS  difference  plots. 

In  general,  the  performance  of  the  Advanced  CIRIS  filter  is  better  than  that  of  the 
current  CIRIS  II  filter.  The  exception  is  in  the  range-rate  measurement  processing.  The 
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Advanced  CIRIS  filter  range-rate  measurement  residuals  are  unacceptably  luge  during 
aircraft  maneuvers,  resulting  in  an  excessive  percentage  of  range-rate  measurements  being 
rejected  by  the  filter. 


The  third  goal  was  the  development  of  a  fixed-interval  smoothing  procedure  compat¬ 
ible  with  the  MSOFE-based  Advanced  CIRIS  filter.  The  description  of  such  a  procedure 
is  given  in  Chapter  VI,  but  it  was  not  implemented  and  tested  due  to  a  shortage  of  time  to 
work  on  it.  The  implementation  of  this  procedure  will  be  significantly  more  difficult  than 
simply  proposing  it.  The  primary  obstacles  are  the  luge  amounts  of  disk  storage  required 
and  the  need  to  ensure  compatability  with  the  state  switching  routine  implemented  in  the 
Advanced  CIRIS  filter. 

7. 2  Recommendations 

The  following  recommendations  ue  divided  into  three  categories:  softwue  develop¬ 
ment,  the  need  for  additional  data,  and  directions  for  future  reseuch. 

7.2.1  Software  Development  The  reseuch  described  in  this  thesis  required  the  de¬ 
velopment  of  a  significant  amount  of  softwue.  In  addition  to  the  MSOFE  user-subroutines 
which  implement  the  Advanced  CIRIS  filter,  there  ue  several  small  utility  programs  needed 
for  data  file  manipulation,  MATRIXx  programs  for  making  plots,  and  operating  system 
command  files  for  running  batch  jobs  and  recompiling  the  filter  subroutines.  Yet,  the  de¬ 
velopment  of  softwue  is  not  the  main  focus  of  the  reseuch.  The  main  focus  is  the  function 
and  performance  of  the  software:  the  Kalman  filter.  However,  it  is  the  author’s  firm  opin¬ 
ion  that  softwue  functionality  and  performance  cannot  be  truly  achieved  in  the  absence 
of  thoughtful  design,  implementation,  and  maintenance  of  the  softwue.  The  recommen¬ 
dation  is  that  future  student  reseuchers  be  encouraged  and  given  credit  for  devoting  a 
larger  portion  of  their  efforts  to  developing  softwue  that  is  robust,  documented,  and  easily 
modified. 

7.2.2  Data  Problems  The  use  of  empirical  data  in  this  reseuch,  though  necessary  to 
validate  the  correctness  of  the  filter  model,  presented  many  problems  and  required  the  use 
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of  certain  approximations.  The  data  tapes  (CIRIS,  laser,  track)  contained  a  small  number 
of  data  dropouts,  duplicate  records,  and  obviously  erroneous  data  values.  The  presence  of 
these  occasional  bad  records  cannot  always  be  ignored;  they  must  be  located  and  either 
“patched”  or  skipped.  Given  the  huge  size  of  the  data  files,  this  is  a  time-consuming  task. 

Even  more  troublesome  than  bad  data  records,  certain  information  required  by  the 
analytic  description  of  the  filter  model  is  simply  unavailable  on  the  data  tapes.  The 
three  items  in  this  category  are  the  raw  INS  accelerometer  outputs,  barometric  altitude 
rate,  and  aircraft  attitude  rates.  The  accelerometer  outputs  are  available  from  the  INS 
computer  and  attempts  have  been  made  to  record  them,  but  the  great  majority  of  the 
recorded  acceleration  values  are  zero  [11].  This  is  possibly  due  to  a  timing  problem  in 
the  recording  process.  In  the  absence  of  this  information,  platform  frame  accelerations 
must  be  “reverse  engineered”  (synthesized)  from  the  recorded  INS  velocities  and  other 
recorded  information.  The  result  is  acceleration  data  that  is  representative  of  the  true 
average  acceleration  between  sampling  times.  Since  this  acceleration  information  is  used 
only  in  defining  the  INS  error  state  dynamics  matrix  for  propagation  of  the  INS  error  state 
estimates,  and  is  not  a  direct  measurement  to  the  filter,  the  use  of  consistent,  synthesized 
accelerations  should  have  a  negligible  impact  on  filter  performance.  The  recommendation 
is  to  investigate  why  the  raw  accelerometer  outputs  cannot  be  reliably  recorded  and  correct 
the  problem. 

The  barometric  altitude  rate  (h\,)  is  needed  only  for  filter  emulation  of  the  internal 
INS  calculation  of  the  vertical  channel  aiding  gains,  as  described  in  Chapter  III.  However, 
this  rate  data  is  not  recorded  (it  may  not  be  directly  available  from  the  CADC)  so  it  must 
be  synthesized  by  differencing  consecutive  barometric  altitude  measurements  and  dividing 
by  the  measurement  interval.  This  explains  the  appearance  of  the  baro-altitude  rate  plot 
in  Figure  B.6.  The  recommendation  is  to  continue  this  procedure  at  a  higher  sampling 
rate  (ideally,  the  same  sampling  rate  used  by  the  INS)  or,  if  available,  directly  record  the 
aiding  gain  values  as  calculated  by  the  INS. 

The  aircraft  attitude  rates  (ip,  9 ,  < p)  are  needed  to  calculate  the  transformation  ma¬ 
trices  for  use  in  the  lever-arm  compensation  procedures  for  the  range-rate  measurements. 
These  three  quantities  are  not  available  on  the  data  tape,  so  they  must  be  synthesized  by 


7-5 


differencing  consecutive  attitude  angles  and  dividing  by  the  sampling  interval.  The  rec¬ 
ommendation  is  to  continue  this  procedure  at  a  higher  sampling  rate  or  to  directly  record 
the  attitude  rates  if  they  are  available  from  the  INS. 

7 .2.3  Future  Research  For  reasons  discussed  in  the  preface,  the  scope  of  this  re¬ 
search  overlapped  that  of  Solomon  [13]  to  a  great  extent.  Differences  include  the  analytic 
development  of  the  lever  arm  compensation  procedure,  the  analytic  description  of  the 
filter  measurement  model  using  extended  Kalman  filter  notation,  modification  of  the  at¬ 
mospheric  propagation  delay  error  model,  and  the  use  of  two  additional  real  data  sets, 
with  independent  position/velocity  measurements,  to  evaluate  filter  performance.  As  a 
result  of  this  additional  work,  the  concept  of  the  Advanced  CIRIS  filter  has  been  more 
completely  characterized  and  tested,  and  so  it  has  greater  potential  for  “operational”  use, 
should  CIGTF  desire  to  use  it.  The  filter  software  is  now  structured  well  enough  to  be  a 
good  starting  point  for  continued  research  at  both  CIGTF  and  AFIT. 

Qualitatively,  the  performance  of  the  Advanced  CIRIS  filter  is  sufficiently  greater 
than  that  of  the  CIRIS  II  filter  to  warrant  further  development.  This  improvement  is 
primarily  due  to  the  filter  structure  (additional  observable  states,  more  accurate  propaga¬ 
tion  and  update  algorithms)  since  the  filter  initial  covariance  values  and  atmospheric  error 
model  parameters  are  not  yet  particularly  well  tuned.  The  effort  to  tune  filter  parameters 
precisely  and  quantify  filter  estimation  accuracy  should  be  continued  at  CIGTF  instead  of 
at  AFIT  because  of  the  availability  of  additional  data  sets  and  other  necessary  resources. 
As  an  aid  in  this  tuning  process,  it  would  be  useful  to  collect  additional  CIRIS  data  from 
tests  designed  to  achieve  high  observability  of  specific  error  states.  The  key  to  correctly 
modeling  the  atmospheric  delay  errors  is  availability  of  power  spectral  density  (and  cross 
spectral  density)  information.  Such  information  is  most  easily  collected  from  a  long  dura¬ 
tion  static  test  involving  severed  carefully  placed  transponders  sampled  at  a  high  (0.1 — 1.0 
Hz),  uniform  sampling  rate. 

Further  research  at  AFIT  should  focus  on  two  areas:  use  of  Global  Positioning  Sys¬ 
tem  (GPS)  data  and  implementation  of  a  fixed-interval  smoother.  The  incorporation  of 
GPS,  particularly  differential  GPS,  measurements  into  the  Advanced  CIRIS  filter  should 
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be  pursued  first.  The  development  of  a  complete,  correct  error  model  for  differential  GPS 
measurements  promises  to  be  the  single  most  effective  way  to  increase  the  estimation  ac¬ 
curacy  of  the  filter.  The  Advanced  CIRIS/GPS  filter  should  have  the  capacity  to  process 
GPS  range  and  range-rate  measurements  from  at  least  four  different  satellites,  but  prefer¬ 
ably  from  all  satellites  whose  measurements  Eire  recorded.  One  of  the  great  advEintages 
of  an  MSOFE  based  filter  is  the  ability  to  quickly  reconfigure  the  filter  dimension  (state 
vector  length)  to  provide  exactly  the  right  number  of  error  states  required  for  the  number 
of  transponders  or  satellites  used  in  a  particular  data  set.1  The  development  of  the  Ad¬ 
vanced  CIRIS/GPS  filter  should  conclude  with  a  performance  evEduation  using  real  data. 
However,  much  of  the  initial  development  could  be  accomplished  using  MSOFE’s  Monte 
Carlo  simulation  mode.  Initisd  development  using  Monte  Carlo  simulation  would  avoid  the 
previously  mentioned  pitfalls  associated  with  empirical  data. 

The  implementation  of  a  fixed-interval  smoother,  either  the  one  described  in  Chap¬ 
ter  VI  or  one  based  on  euti  extensively  modified  version  of  MSOFE,  should  be  undertsiken 
to  determine  the  potentiEil  increase  in  estimation  accuracy  over  the  stEind-alone,  forward- 
pass  filter.  In  either  cEise,  initial  development  should  use  Monte  Carlo  simulation  to  permit 
a  precise,  quantitative  assessment  on  the  accuracy  of  the  smoothed  estimates  versus  the 
forward-pass  estimates.  The  development  of  a  fixed-interval  smoothing  capability  for 
MSOFE  would  be  eui  extremely  useful  achievement,  but  this  is  a  topic  worthy  of  an  inde¬ 
pendent  thesis  effort.  If  attempted,  it  should  be  done  in  close  cooperation  with  the  authors 
of  MSOFE. 


‘This  is  the  optimal  alternative  to  the  ad  hoc  state  sharing  procedure  currently  used  with  the  transpon¬ 
ders,  but  it  may  result  in  a  large  state  vector  length. 
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Appendix  A.  Initial  Covariance  Values  and  Filter  Matrix  Definitions 


A-l 


Table  A.l.  Non-Zero  Elements  of  Pq 


Element 

State 

Definition 

Value  ( lcr) 

1,1 

SQ  x 

Magnitude  of  error  between  the  true  frame’s  X 
axis  and  the  computer  frame’s  X  axis. 

50  arcsecs 

2,2 

8Qy 

Magnitude  of  error  between  the  true  frame’s  Y 
axis  and  the  computer  frame’s  Y  axis. 

50  arcsecs 

3,3 

6®z 

Magnitude  of  error  between  the  true  frame’s  Z 
axis  and  the  computer  frame’s  Z  axis. 

50  arcsecs 

4,4 

4>x 

Magnitude  of  error  between  the  true  frame’s  X 
axis  and  the  platform  frame’s  X  axis. 

100  arcsecs 

5,5 

4>y 

Magnitude  of  error  between  the  true  frame’s  Y 
axis  and  the  platform  frame’s  Y  axis. 

100  arcsecs 

6,6 

4>z 

Magnitude  of  error  between  the  true  frame’s  Z 
axis  and  the  platform  frame’s  Z  axis. 

180  arcmins 

8VX 

Magnitude  of  error  between  the  true  frame  X 
velocity  and  the  computer  frame  X  velocity. 

1.0  ft/sec 

8,8 

SVy 

Magnitude  of  error  between  the  true  frame  Y 
velocity  and  the  computer  frame  Y  velocity. 

1.0  ft/sec 

9,9 

6VZ 

Magnitude  of  error  between  the  true  frame  Z 
velocity  and  the  computer  frame  Z  velocity. 

1.0  ft/ sec 

10,10 

8h 

Magnitude  of  error  between  the  true  altitude 
and  the  computer  altitude. 

100  ft 

m 

Magnitude  of  error  between  the  true  altitude 
and  the  computer  altitude  (1  sec  delay). 

100  ft 

12,12 

SS3 

Magnitude  of  vertical  channel  aiding  state. 

1  ft 

13,13 

SS4 

Magnitude  of  vertical  channel  aiding  state. 

All  other  elements  are  zero. 
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Table  A. 2.  Non-Zero  Elements  of  P0,  Continued 


Element 

State 

Definition 

Value  (lcr) 

14,14 

Bxc 

Magnitude  of  gyro  correlated  drift  rate  along  X 
axis. 

0.002  deg / hr 

15,15 

Byc 

16,16 

Bzc 

Magnitude  of  gyro  correlated  drift  rate  along  Z 
axis. 

0.005  deg/hr 

17,17 

V.Y  c 

Magnitude  of  accelerometer  and.  velocity  quan¬ 
tizer  correlated  noise  along  X  axis. 

3  jj.g 

18,18 

Vyc 

Magnitude  of  accelerometer  and  velocity  quan¬ 
tizer  correlated  noise  along  Y  axis. 

3  ng 

19,19 

v*c 

Magnitude  of  accelerometer  and  velocity  quan¬ 
tizer  correlated  noise  along  Z  axis. 

3  ng 

20,20 

Magnitude  of  gravity  vector  error  along  X  axis. 

35  pg 

21,21 

6Gy 

Magnitude  of  gravity  vector  error  along  Y  axis. 

35  ng 

22,22 

6GZ 

Magnitude  of  gravity  vector  error  along  Z  axis. 

35  /j.g 

23,23 

Bx 

Magnitude  of  gyro  drift  rate  bias  repeatability 
along  X  axis. 

24,24 

By 

Magnitude  of  gyro  drift  rate  bias  repeatability 
along  Y  axis. 

0.0045  deg/hr 

25,25 

Bz 

Magnitude  of  gyro  drift  rate  bias  repeatability 
along  Z  axis. 

0.01  deg/hr 

All  other  elements  are  zero. 
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Table  A. 3.  Non-Zero  Elements  of  Pq,  Continued 


mm 

Definition 

Value  (l<r)  | 

26,26 

Magnitude  of  baro-altimeter  correlated  bias 
noise. 

100  ft 

27,27 

f>hb2 

Magnitude  of  baro-altimeter  constant  bias  er¬ 
ror. 

28,28 

_ 

Magnitude  of  baro-altimeter  scale  factor  error. 

_ 

4  % 

Ail  other  elements  are  zero. 
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Table  A. 4.  Non-Zero  Elements  of  Pq,  Continued 


Element  i  State 


Definition 


29,29 

5Erc 

Magnitude  of  the  transponder  interogator  range 
measurement  calibration  error. 

1  ft 

30,30 

8  Erc 

— 

Magnitude  of  the  transponder  interogator 
range-rate  measurement  calibration  error. 

0.01  ft/se  c 

I  31,31  8 Magnitude  of  transponder  1  position  error  5  ft 

j  ( XE  axis). 


32,32 

8YTl 

— 

Magnitude  of  transponder  1  position  error 
(YE  axis). 

— 

5  ft 

33,33  8Zpt  Magnitude  of  transponder  1  position  error  5  ft 
( ZE  axis). 


34,34  8Ati  Magnitude  of  transponder  1  atmospheric  scale 
factor  error  (line-of-sight). 


35,35  8Xt2  Magnitude  of  transponder  2  position  error  5  ft 
(XE  axis). 


Magnitude  of  transponder  2  position  error 
(YE  axis). 


37,37  8Zt2  Magnitude  of  transponder  2  position  error  5  ft 
)  (ZE  axis). 


38,38  8At3  Magnitude  of  transponder  2  atmospheric  scale  10  ppm 
factor  error  (line-of-sight). 


39,39 


40,40 


41,41 


42,42 


8Yp3  Magnitude  of  transponder  3  position  error 
axis). 


6Zj3  Magnitude  of  transponder  3  position  error 
( ZE  axis). 


8At3  Magnitude  of  transponder  3  atmospheric  scale 
factor  error  (line-of-sight). 


Table  A. 5.  Non-Zero  Elements  of  Po,  Continued 


Definition 

Value  (lcr) 

43,43 

SX  r4 

Magnitude  of  transponder  4  position  error 
(XE  axis). 

5  ft 

44,44 

6YT< 

Magnitude  of  transponder  4  position  error 
( Ye  axis). 

5  ft 

45,45 

SZj4 

Magnitude  of  transponder  4  position  error 
( ZE  axis). 

5  ft 

46,46 

Magnitude  of  transponder  4  atmospheric  scale 
factor  error  (line-of-sight). 

10  ppm 

47,47 

SXTb 

Magnitude  of  transponder  5  position  error 
(XE  axis). 

5  ft 

48,48 

6Yt> 

Magnitude  of  transponder  5  position  error 
(Ye  axis). 

5  ft 

49,49 

Magnitude  of  transponder  5  position  error 
( ZE  axis). 

5  ft 

50,50 

6  At,, 

Magnitude  of  transponder  5  atmospheric  scale 
factor  error  (line-of-sight). 

10  ppm 

51,51 

6Xn 

Magnitude  of  transponder  6  position  error 
( XE  axis). 

5  ft 

52,52 

SYt6 

Magnitude  of  transponder  6  position  error 
(Ye  axis). 

•5  ft 

53,53 

6ZT<i 

Magnitude  of  transponder  6  position  error 
( ZE  axis). 

5  ft 

54,54 

sat6 

Magnitude  of  transponder  6  atmospheric  scale 
factor  error  (line-of-sight). 

10  ppm 

Magnitude  of  transponder  7  position  error 
( XE  axis). 

5  ft 

Magnitude  of  transponder  7  position  error 
(Ye  axis). 

5  ft 

Magnitude  of  transponder  7  position  error 
( ZE  axis). 

5  ft 

58,58 

8At7 

Magnitude  of  transponder  7  atmospheric  scale 
factor  error  (line-of-sight). 

10  ppm 

Table  A. 6.  Non-Zero  Elements  of  Po,  Continued 


Element 

59,59 


60,60 


61,61 


62,62 


State 


6YTt 


SZTi 


6ATft 


Definition 

Magnitude  of  transponder  8  position  error 
( XE  axis). 

Magnitude  of  transponder  8  position  error 
( Ye  axis). 

Magnitude  of  transponder  8  position  error 
( ZE  axis). 

Magnitude  of  transponder  8  atmospheric  scale 
factor  error  (line-of-sight). 


63,63 


64,64 


65,65 


66,66 


67,67 


68,68 


6XT, 


6Yt, 


6At« 


SXTv 


syTii 


Magnitude  of  transponder  9  position  error 
(XE  axis). 

Magnitude  of  transponder  9  position  error 
(Ye  axis). 

Magnitude  of  transponder  9  position  error 
( ZE  axis). 

Magnitude  of  transponder  9  atmospheric  scale 
factor  error  (line-of-sight). 

Magnitude  of  transponder  10  position  error 
( XE  axis). 

Magnitude  of  transponder  10  position  error 
(Ye  axis). 


69,69 


70,70 


6Zt10  Magnitude  of  transponder  10  position  error 
( ZE  axis). 

tfAru,  Magnitude  of  transponder  10  atmospheric  scale 
factor  error  (line-of-sight). 


All  other  elements  are  zero. 


Value  ( 1  cr  ) 

Vft 


5  ft 


5  ft 


10  ppm 


5  ft 


5  ft 


5  ft 


10  ppm 


5  ft 


5  ft 


5  ft 


10  ppm 


Table  A. 7.  Non-Zero  Elements  of  F 


Element  I  Symbol  or  Formula  |  Value 


—u>EPY 


-Cry 


UEPX 


-UJEPX 


U>IEY 


UIPZ 


UIEZ 


-UtE. X 


Uipz 


-Uipy 


Crx 


1.0 


1.0 


-^lEY 


VIEX 


U>ipy 


-UJjpX 


1.0 


1.0 


-2(Vy  &IEY  +  Vz^iez) 


2{Vy<j)IEx) 


2{Vzuiex) 


-A  z 


Ay 


c 

RX 

t 

UJEPY 

t 

7,8 

2  &IEZ 

t 

7,9 

-{UEPY  +  2  ujiey) 

t 

7,17 

1.0 

1.0 

7,20 

1.0 

1.0 

and/or  trajectory  dependent 


Table  A. 8.  Non-Zero  Elements  of  F,  Continued 


Element 

Symbol  or  Formula 

8,1 

2  Vx  uiey 

8,2 

-2(  Vx  u>rEX  +  Vz  uiez) 

8,3 

2  Vz  wjey 

8,4 

Az 

8,6 

-AX 

8,7 

“2  uiez 

Value 


-Vz  Cry 


UEPX  +  2  U//EX 


2  Vx  uiez 


2  Vy  uiez 


-2(  Vy  ujiey  +  Vx  UIEX) 


- Ay 


UEPY  +  2  UJIEY  +  Vx  Crx 


9,8 

-UIEPX  -  2  >^IEY  +  Vy  Cry 

9,10 

2  Go/ A 

9,11 

-k2 

9,12 

-1.0 

9,13 

k2 

9,19 

1.0 

9,22 

1.0 

9,23 

k2 

9,27 

k2 

9,28 

h  K2 

10,9 


10,11 


10,13 


10,23 


10,27 


10,28 


Ki  -  1.0 


K  i 


K  x 


hKx 


11,10 

1.0 

bxh 

11,11 

-1.0 

-1.0 

12,11 

Kz 

t 

12,13 

-k3 

t 

12,23 

-k3 

t 

12,27 

-k3 

t 

12,28 


time  and/or  trajectory  dependent 
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Table  A. 9.  Non-Zero  Elements  of  F,  Continued 


Element 


13,10 


13,11 


13,13 


13,23 


13,28 


14,14 


15,15 


16,16 


17,17 


18,18 


19,19 


20,20 


21,21 


22,22 


23,23 


34,34 


38,38 


42,42 


46,46 


50,50 


54,54 


58,58 


62,62 


66,66 


70,70 


Symbol  or  Formula 


K, 


-Ka 


A  4  —  1.0 


A*4  fishr. 


-KAVZ 


~Pb 


-Pb, 


as- 


-pBZr 


-Pv 


As- 


-Pv-r, 


-Pv 


Zl±- 


p6Gy 


P$  G  y 


”  p6  G  2 


-Pshr 


-1/Tal 


-1/t02 


-1/Ta3 


-l/fa4 


-1/Ta5 


-1/Ta6 


-l/T-gT 


-1/Ta8 


-1/Ta9 


-1/7-alO 


f  time  and/or  trajectory  dependent 


Value 


t 


t 


-1/300  (s-1) 


-1/300  (s-1) 


-1/300  (s'1) 


-1/600  (s"1) 


-1/600  (s'1) 


-1/600  (s-1) 


f 


t 


-1/600  (s 


-l  ’ 


-1/300  (s'1) 


-1/300  (s"1) 


-1/300  (s'1) 


-1/300  (s'1) 


-1/300  (s 


-l  i 


-1/300  (s'1) 


-1/300  (s-1) 


-1/300  (s31") 


-1/300  (s'1) 


-1/300  (s'1) 


Table  A. 10.  Non-Zero  Elements  of  Q 


Element  Symbol  or  Formula 


Value 


9.3  -  10~9  (r/s3) 


9.3-10  ~9(f2/s3) 


2  •  10_lo/300  (a-1) 


2  •  10_1°/300  (a"1) 


2  •  10-lo/300  (a"1) 


2  •  10~lo/300  (s'1) 


2  •  10"lo/300  (s_1) 


2  •  10~lo/300  (s'1) 


66,66 

2  <r*9/ra9 

1  2  •  10_lo/300  (s'1)  1 

70,70 

2o-aio/raio 

t  time  and/or  trajectory  dependent 


The  off-diagonal  elements  indexed  by  all  possible 
pairs  of  the  form  (i,j),  i  ^  j,  selected  from  {34, 
38,  42,  46,  50,  54,  58,  62,  66,  70}  are  calculated  as 


Qi.j  = 


2  P  &ai  &a 


'Ta{  Taj 


where  p  is  the  correlation  coefficient.  The  value 
used  for  p  is  0.8. 


Table  A. 11.  Non-Zero  Elements  of  Hr 


Element 

U 

S&X 

Hi 

1,2 

6®Y 

Hi 

1,10 

Sh 

Hi 

1,29 

6  Ere 

jBjfPWl 

h#XT] 

sxT] 

Xt,~Xa 

R, 

1  ,#Yji 

6YT, 

?T  -Ya 
j _ ii _ ± 

' 

1  #zT, 

6Zt, 

Zx  ~~Z  A 

H — 

R, 

l.iMr, 

SAt, 

r3 

All  other  elements  are  zero. 

NOTE  Is  The  H'  transformation  projects  the  INS  position  errors,  as  defined  in  the  state 
vector,  into  the  E-frame.  This  transformation  is  required  because  the  range  measurements 
are  defined  in  terms  of  E-frame  coordinates.  The  elements  of  H'  are  calculated  as 


Hi  Hi  Hi 


xT, -*a 

- Sc: - 

Rj 


where  Cjy  is  defined  in  Equation  (2.10)  and  T% 


?t,~Ya  Zt,~Za 

R,  R, 


C%T 


N 

X 


is  defined  in  Equation  (3.5). 


NOTE  2:  The  #  symbol  represents  the  number  of  the  corresponding  state  in  the  10 
transponder  window  for  the  transponder  being  sampled  at  the  applicable  measurement 
update  time. 
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Table  A.  12.  Non- Zero  Elements  of  Hr 


Element 

State 

Formula 

1,1 

6®x 

Hi' 

1,2 

6Qy 

H'2' 

1,10 

6h 

Hi' 

1,7 

SVX 

Hi' 

1,8 

SVY 

Hi' 

1,9 

SVZ 

Hi' 

1,30 

SErc 

1.0 

All  other  elements  are  zero. 

NOTE  Is  The  H"  transformation  projects  the  INS  position  and  velocity  errors,  as  defined 
in  the  state  vector,  into  the  E-frame.  This  transformation  is  required  because  the  range- 
rate  measurements  are  defined  in  terms  of  E-frame  coordinates.  The  elements  of  H"  are 
calculated  as 

H"  Yf"  LT"  = 

1  tl2  ii3 


(At,  -  XA)-XA  +  (Xtj  -  XA)(?Tj  -  Ya)Ya  +  (AT,  -  XA)(ZTi  -  ZA)ZA  -  R$. XA 


4 


(?Tt  -  Ya  )(At,  -  XA)XA  +  (fT,  -  YA)2YA  +  (YTj  -  ?A)(ZT,  -  ZA)ZA  -  RfYA 

4 

(ZT,  -  ZA)(XTj  -  A*)!*  +  ( ZTj  -  Za)(YTj  -  YA)YA  +  (ZTj  -  ZA?ZA  - 


4 


E 


C%T 


N  *  X 


and 


XI  ft  u"  U" 

n4  ii5  xl6 

= 

,Tt  -X a  Yt-Ya  Z? -ZA 

— b*. - - bz -  —  —  - 

H,  R,  Rj 


c 


E 

P 


where  is  defined  in  Equation  (2.10),  T %  is  defined  in  Equation  (3.5),  and  Cp  is  defined 
in  Equation  (2.11) 
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Table  A.  13.  Non- Zero  Elements  of 


Element 

1,10 

Sh 

-1.0 

1,26 

Sh(,x 

1.0 

1,27 

Shb2 

1.0  j 

1,28 

Shf,  3 

h, 

All  other  elements  are  zero. 
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Table  A.  14.  Non- Zero  Elements  of  R 


Element 

Symbol  or  Formula 

Value 

u 

<4 

■l.lltiW 

2,2 

0.25  {f2/s2) 

3,3 

< 

9.0  (f2) 

In  the  MSOFE-based  filter,  R  never  appears  as 
a  3-by-3  matrix.  MSOFE’s  scalar  update  pro¬ 
cedure  sequentially  uses  the  diagonal  elements  of 
the  3-by-3  R  as  three  separate  1-by-l  R  ma¬ 
trices. 
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Appendix  B.  Recorded  Data  for  Track  Run  E 


B-l 


54800 


55000  55200  55400  55600  55800 

IRIG  TIME  (seconds) 


Figure  B.l.  TDAS  Position,  Track  Run  E 


Figure  B.2.  TDAS  Velocity,  Track  Run  E 
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Figure  B.4.  INS  Indicated  Velocity,  Track  Run  E 
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Figure  B.5.  INS  Indicated  Acceleration  and  Wander  Angle,  Track  Run  E 
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Figure  B.6.  Baro-altimeter  Altitude  and  Altitude  Rate,  Track  Run  E 
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Figure  B.7.  Transponder  005  Range  and  Range-Rate,  Track  Run  E 
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Figure  B.8.  Transponder  102  Range  and  Range-Rate,  Track  Run  E 
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Figure  B.9.  Transponder  181  Range  and  Range-Rate,  Track  Run  E 
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Figure  B.10.  Transponder  211  Range  and  Range-Rate,  Track  Run  E 
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Figure  B.ll.  Transponder  212  Range  and  Range-Rate,  Track  Run  E 
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Figure  B.12.  Transponder  216  Range  and  Range-Rate,  Track  Run  E 
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Appendix  C.  Recorded  Trajectory  Data,  Yuma  Flight 
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Figure  C.l.  Laser  Ranging  System  Indicated  Position,  Yuma  Flight 
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Figure  C.2.  Laser  Ranging  System  Indicated  Velocity,  Yuma  Flight 
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Figure  C.4.  LN-39  INS  Indicated  Velocity,  Yuma  Flight 
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Figure  C.5.  INS  Indicated  Acceleration  and  Wander  Angle,  Yuma  Flight 
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Abstract 

The  Completely  Integrated  Reference  Instrumentation  System  (CIRIS)  is  operated 
bv  the  Central  Inertial  Guidance  Test  Facility  (CIGTF),  at  Holloman  AFB,  New  Mex¬ 
ico.  CIRIS  functions  as  a  reference  navigation  system  used  for  evaluating  the  accuracy 
of  other  types  of  navigations  systems.  As  a  reference  standard,  the  root-mean-square 
errors  in  CIRIS  estimates  of  aircraft  trajectory  variables  must  be  maintained  an  order  of 
magnitude  smaller  than  the  corresponding  errors  of  the  navigation  system  under  test.  The 
primary  hardware  components  of  CIRIS  are  a  reference  inertial  navigation  system  (INS), 
a  baro-altimeter,  an  array  of  ground-based  transponders,  a  transponder  interrogator,  and 
data  recording  equipment.  The  transponder  equipment  provides  transponder-to-aircraft 
range  and  range-rate  measurements  during  test  flights.  The  primary  software  compo¬ 
nent  of  CIRIS  is  a  Kalman  filter  program  which  processes  the  recorded  measurements  and 
estimates  the  true  position  and  velocity  of  the  aircraft  throughout  the  test  flight. 

The  estimation  accuracy  of  CIRIS  must  be  increased  so  that  CIRIS  can  serve  as  a 
benchmark  for  measuring  the  accuracy  of  Global  Positioning  System  (GPS)  aided  inertial 
navigation  systems.  This  thesis  documents  the  continuation  of  research  to  develop  com¬ 
pletely  new  Kalman  filter  software  for  CIRIS.  A  70-state  filter,  based  on  the  Multimode 
Simulation  for  Optimal  Filter  Evaluation  (MSOFE)  program,  developed  in  previous  re¬ 
search  is  the  starting  point.  This  70-state  filter  models  error  dynamics  associated  with  the 
CIRIS  Litton  LN-39  INS,  baro-altimeter,  and  transponder  equipment.  In  this  thesis,  the 
model  for  atmospheric  effects  on  transponder  range  measurements  is  refined  and  the  filter 
is  modified  to  process  barometric  altitude  measurements  in  addition  to  the  transponder 
measurements.  The  performance  of  the  resulting  filter  is  evaluated  using  real  CIRIS  data 
recorded  during  a  slow  speed  ground  test  and  an  aircraft  flight  test.  The  filter  position  and 
velocity  estimates  are  compared  to  independent  measurements  of  the  same  quantities.  The 
structure  for  a  companion  fixed-interval  smoother  program  is  proposed  and  discussed,  hut 
not  implemented.  Future  research  is  expected  to  extend  the  filter  to  process  differential 
mode  GPS  measurements. 
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